Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Alternatives/improvements to mldivide?

Subject: Alternatives/improvements to mldivide?

From: PLH

Date: 8 Nov, 2010 16:25:20

Message: 1 of 10

I'm trying to speed up my code, and the current bottleneck is in using the BACKSLASH operator (mldivide).

Essentially, the relevant part of my code looks like:

N = 6000000; % a big number
for j=1:N
    y_new = A\y_old
    A = my_func(y_old); % uses an optimised (much reduced) version of spdiags

The variables y_new and y_old are Jx1 vectors, and the matrix A is square (JxJ) and banded (the 5 inner diagonals), but with three nonzero elements in each of the opposite corners. By this I mean,

A(1,J)
A(1,J-1)
A(2,J)
A(J,1)
A(J,2)
A(J-1,1)

are nonzero. So, A is almost pentadiagonal, but really it has 9 bands (5 in the middle, 2 at the top, 2 at the bottom). Put another way, the linear system I'm solving is pentadiagonal with periodic boundary conditions imposed on the solution y_new.

I store A as a sparse matrix, and use mldivide. It works, but I'm looking for speedups. I've tried "factorize", which is offered on the file exchange, but this almost 1.5 times slower.

If someone could suggest how to obtain a speedup, I would be extremely grateful. Whether that's a straight alternative to mldivide, or whether someone out there knows an algorithm for the particular class of matrix I have, I'm open to all suggestions.

Thanks,

PLH.

Subject: Alternatives/improvements to mldivide?

From: Steven_Lord

Date: 8 Nov, 2010 19:11:47

Message: 2 of 10



"PLH " <paulhalkyard@googlemail.com> wrote in message
news:ib989g$fmr$1@fred.mathworks.com...
> I'm trying to speed up my code, and the current bottleneck is in using the
> BACKSLASH operator (mldivide).
>
> Essentially, the relevant part of my code looks like:
>
> N = 6000000; % a big number
> for j=1:N
> y_new = A\y_old
> A = my_func(y_old); % uses an optimised (much reduced) version of
> spdiags

Are you CERTAIN that the bottleneck is \ instead of my_func?

It looks almost like you're trying to build your own iterative solver.

> The variables y_new and y_old are Jx1 vectors, and the matrix A is square
> (JxJ) and banded (the 5 inner diagonals), but with three nonzero elements
> in each of the opposite corners. By this I mean,
> A(1,J)
> A(1,J-1) A(2,J)
> A(J,1)
> A(J,2)
> A(J-1,1)
>
> are nonzero. So, A is almost pentadiagonal, but really it has 9 bands (5
> in the middle, 2 at the top, 2 at the bottom). Put another way, the
> linear system I'm solving is pentadiagonal with periodic boundary
> conditions imposed on the solution y_new.
>
> I store A as a sparse matrix, and use mldivide. It works, but I'm looking
> for speedups. I've tried "factorize", which is offered on the file
> exchange, but this almost 1.5 times slower.

How long does it currently take, and how long would you like it to take?
[I'm trying to see whether you're hoping for a level of performance that is
feasible to achieve or if you're looking for, essentially, magic.]

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: Alternatives/improvements to mldivide?

From: PLH

Date: 9 Nov, 2010 09:02:03

Message: 3 of 10

Hello and thanks for the reply.

I identified the bottleneck using Matlab's profiler, which states that 60% of the run-time is spent in mldivide, 15% is spent in "sparse", and the other 25% is "all other lines" (integrations and the like). Mostly, I've replaced a few Matlab functions (such as spdiags, hence the usage of sparse) with very specific equivalents that are written inline to achieve some speedup. Without using a new approach, I believe the only unit of my code left to attack is mldivide, but I'd appreciate any suggestions along lines too!

Currently, each loops runs for about 10 minutes, which in itself is fine. However, I want to repeat the process for many different parameter sets (the matrix A is a function of 5 parameters set at run-time), so any small speedup will be noticeable - I'm not expecting order of magnitude changes, but a factor of 2 would be nice.

Yes, you're correct - it's an iterativer solver (Crank-Nicolson with some bells and whistles for a nonlinear PDE). If nothing can be done, then it's not critical - if it was, I'd just spend the next 20 coffees writing it in Fortran. But if someone out there knew the XYZ algorithm for just such a class of banded matrix, or knew some way of informing mldivide of what algorithm to use, then it would have been useful to know.

Any suggestions, or should I be happy with what I have?

Subject: Alternatives/improvements to mldivide?

From: Steven_Lord

Date: 9 Nov, 2010 15:04:48

Message: 4 of 10



"PLH " <paulhalkyard@googlemail.com> wrote in message
news:ibb2mb$4s4$1@fred.mathworks.com...
> Hello and thanks for the reply.
>
> I identified the bottleneck using Matlab's profiler, which states that 60%
> of the run-time is spent in mldivide, 15% is spent in "sparse", and the
> other 25% is "all other lines" (integrations and the like). Mostly, I've
> replaced a few Matlab functions (such as spdiags, hence the usage of
> sparse) with very specific equivalents that are written inline to achieve
> some speedup. Without using a new approach, I believe the only unit of my
> code left to attack is mldivide, but I'd appreciate any suggestions along
> lines too!

Well, you only showed a small segment of your code, so it would be difficult
to offer suggestions for other avenues of attack.

> Currently, each loops runs for about 10 minutes, which in itself is fine.
> However, I want to repeat the process for many different parameter sets
> (the matrix A is a function of 5 parameters set at run-time), so any small
> speedup will be noticeable - I'm not expecting order of magnitude changes,
> but a factor of 2 would be nice.

Six million iterations through the loop (based on your posted value of N) in
600 seconds is ten thousand iterations per second or one iteration per one
ten-thousandth of a second. That's pretty quick.

Are you iterating that many times through the loop because that's an upper
bound on how long it'll take to converge or do you really need to iterate
that many times? If you're just looking for convergence, then use a WHILE
loop where one of the conditions it tests is your convergence criterion and
a second condition is the number of iterations (so you don't iterate forever
if the problem diverges or fails to converge for whatever reason.)

while numberOfIterations <= threshold && <some expression testing
convergence>
    % solve
end

> Yes, you're correct - it's an iterativer solver (Crank-Nicolson with some
> bells and whistles for a nonlinear PDE). If nothing can be done, then
> it's not critical - if it was, I'd just spend the next 20 coffees writing
> it in Fortran. But if someone out there knew the XYZ algorithm for just
> such a class of banded matrix, or knew some way of informing mldivide of
> what algorithm to use, then it would have been useful to know.
>
> Any suggestions, or should I be happy with what I have?

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: Alternatives/improvements to mldivide?

From: PLH

Date: 10 Nov, 2010 08:36:03

Message: 5 of 10

I might have been a bit sloppy with my posts so I'll try to clear things up:

Re: small code segment
I only showed a small segment of code because I'm only seeking general advice about

1. The existence of a specific algorithm for solving the class of problem I presented (lots of linear algebra experts on here so I thought that was a reasonable thing to ask).

2. Direct alternatives to mldivide, such as factorize.m (which, unfortunately, is actually slower for my problem).

3. Whether or not it is possible to inform the mldivide function (or easily modify it) to save it having to decide which algorithm to use every iteration.

To be honest, now I think 1 and 2 are a bit optimistic, and instead I'm hoping someone knows something along the lines of 3 (which I only mentioned in subsequent posts).

Re: Iterations
I'm afraid the "6 million" was from memory (since I didn't think being overly specific mattered for my enquiry) - it shouild actually have been more like 1 million. I was just trying to get the point across that I was using mldivide alot, and that it was necessary to rebuild my matrices every time step using the result of the previous iteration.

Re: while loop
Unfortunately, the code isn't aiming to converge to a single function. It's calculating the y(x,t) for prescribed y(x,t=0) -- x is space, t is time. So, I should also have said I need to save each vector every n-th time step (every n-th rather than every 1 in order to save memory, but this takes up less than 3-4% of the runtime according to the profiler so I didn't see it as an issue).

So:
I think I may have made the impression that mldivide's speed was a problem - it's not, but since I'm definitely not a Matlab expert, I presumed I hadn't written my code in the optimum way, so I was systematically trying to improve it. I had already seen massive improvements by replacing spdiags, and the profiler results made it obvious that it was time to try with mldivide - except mldivide is built-in and I can't just open it and start deleting if statements....

I would love to hear any suggestions if there are any, but otherwise I'll just hope that in future releases Mathworks allows mldivide to take some options to tell it what not to try...

Thanks for the help anyway Steve.

PLH

"Steven_Lord" <slord@mathworks.com> wrote in message <ibbnug$2rt$1@fred.mathworks.com>...
>
>
> "PLH " <paulhalkyard@googlemail.com> wrote in message
> news:ibb2mb$4s4$1@fred.mathworks.com...
> > Hello and thanks for the reply.
> >
> > I identified the bottleneck using Matlab's profiler, which states that 60%
> > of the run-time is spent in mldivide, 15% is spent in "sparse", and the
> > other 25% is "all other lines" (integrations and the like). Mostly, I've
> > replaced a few Matlab functions (such as spdiags, hence the usage of
> > sparse) with very specific equivalents that are written inline to achieve
> > some speedup. Without using a new approach, I believe the only unit of my
> > code left to attack is mldivide, but I'd appreciate any suggestions along
> > lines too!
>
> Well, you only showed a small segment of your code, so it would be difficult
> to offer suggestions for other avenues of attack.
>
> > Currently, each loops runs for about 10 minutes, which in itself is fine.
> > However, I want to repeat the process for many different parameter sets
> > (the matrix A is a function of 5 parameters set at run-time), so any small
> > speedup will be noticeable - I'm not expecting order of magnitude changes,
> > but a factor of 2 would be nice.
>
> Six million iterations through the loop (based on your posted value of N) in
> 600 seconds is ten thousand iterations per second or one iteration per one
> ten-thousandth of a second. That's pretty quick.
>
> Are you iterating that many times through the loop because that's an upper
> bound on how long it'll take to converge or do you really need to iterate
> that many times? If you're just looking for convergence, then use a WHILE
> loop where one of the conditions it tests is your convergence criterion and
> a second condition is the number of iterations (so you don't iterate forever
> if the problem diverges or fails to converge for whatever reason.)
>
> while numberOfIterations <= threshold && <some expression testing
> convergence>
> % solve
> end
>
> > Yes, you're correct - it's an iterativer solver (Crank-Nicolson with some
> > bells and whistles for a nonlinear PDE). If nothing can be done, then
> > it's not critical - if it was, I'd just spend the next 20 coffees writing
> > it in Fortran. But if someone out there knew the XYZ algorithm for just
> > such a class of banded matrix, or knew some way of informing mldivide of
> > what algorithm to use, then it would have been useful to know.
> >
> > Any suggestions, or should I be happy with what I have?
>
> --
> Steve Lord
> slord@mathworks.com
> comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
> To contact Technical Support use the Contact Us link on
> http://www.mathworks.com

Subject: Alternatives/improvements to mldivide?

From: PLH

Date: 10 Nov, 2010 09:02:04

Message: 6 of 10

Aha! In case it's useful for anyone else:

http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6TYJ-4Y41V0R-3&_user=2195977&_coverDate=02%2F28%2F2010&_rdoc=1&_fmt=high&_orig=search&_origin=search&_sort=d&_docanchor=&view=c&_searchStrId=1534909065&_rerunOrigin=google&_acct=C000035218&_version=1&_urlVersion=0&_userid=2195977&md5=dd81b2a6736db70e1ace92cc0c850865&searchtype=a

The above is a link to a paper detailing an algorithm designed to obtain the inverse (LU decomposition) for particular class of problem I stated above ("cyclic pentadiagonal" or "periodic pentadiagonal").

PLH

Subject: Alternatives/improvements to mldivide?

From: Steven_Lord

Date: 10 Nov, 2010 15:01:01

Message: 7 of 10



"PLH " <paulhalkyard@googlemail.com> wrote in message
news:ibdlhi$7kb$1@fred.mathworks.com...
> I might have been a bit sloppy with my posts so I'll try to clear things
> up:
>
> Re: small code segment
> I only showed a small segment of code because I'm only seeking general
> advice about
>
> 1. The existence of a specific algorithm for solving the class of problem
> I presented (lots of linear algebra experts on here so I thought that was
> a reasonable thing to ask).
>
> 2. Direct alternatives to mldivide, such as factorize.m (which,
> unfortunately, is actually slower for my problem).
>
> 3. Whether or not it is possible to inform the mldivide function (or
> easily modify it) to save it having to decide which algorithm to use every
> iteration.
>
> To be honest, now I think 1 and 2 are a bit optimistic, and instead I'm
> hoping someone knows something along the lines of 3 (which I only
> mentioned in subsequent posts).

No, it's not possible to do what you suggested. There is a function
LINSOLVE that does allow you to provide some information about the
coefficient matrix you're trying to solve, but that doesn't currently
support sparse matrices.

> Re: Iterations
> I'm afraid the "6 million" was from memory (since I didn't think being
> overly specific mattered for my enquiry) - it shouild actually have been
> more like 1 million. I was just trying to get the point across that I
> was using mldivide alot, and that it was necessary to rebuild my matrices
> every time step using the result of the previous iteration.

Are the changes to the coefficient matrix small enough that the solution
from the previous iteration is "close" to a solution to the modified
problem? If so, perhaps one of the iterative solvers (like GMRES, etc.)
listed in HELP SPARFUN may be useful to you as you can specify the previous
iteration's solution as the starting point for the iterative solver's work
on the next iteration's A.

> Re: while loop
> Unfortunately, the code isn't aiming to converge to a single function.
> It's calculating the y(x,t) for prescribed y(x,t=0) -- x is space, t is
> time. So, I should also have said I need to save each vector every n-th
> time step (every n-th rather than every 1 in order to save memory, but
> this takes up less than 3-4% of the runtime according to the profiler so I
> didn't see it as an issue).
>
> So:
> I think I may have made the impression that mldivide's speed was a
> problem - it's not, but since I'm definitely not a Matlab expert, I
> presumed I hadn't written my code in the optimum way, so I was
> systematically trying to improve it. I had already seen massive
> improvements by replacing spdiags, and the profiler results made it
> obvious that it was time to try with mldivide - except mldivide is
> built-in and I can't just open it and start deleting if statements....
>
> I would love to hear any suggestions if there are any, but otherwise I'll
> just hope that in future releases Mathworks allows mldivide to take some
> options to tell it what not to try...

I'll be honest and say that's unlikely -- it's more likely that LINSOLVE
would be enhanced to support sparse matrices. If you want to "encourage"
that, submit this as an enhancement request via Technical Support so they
can capture your interest in our enhancement database.

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: Alternatives/improvements to mldivide?

From: PLH

Date: 10 Nov, 2010 15:27:03

Message: 8 of 10

I will do just that!

Thanks for your time, I appreciate it!

Subject: Alternatives/improvements to mldivide?

From: Andrew Telatnik

Date: 11 Nov, 2010 14:46:03

Message: 9 of 10

> I store A as a sparse matrix, and use mldivide. It works, but I'm looking for speedups. I've tried "factorize", which is offered on the file exchange, but this almost 1.5 times slower.

I've faced the same performance problems of mldivide, as well as mtimes and other operators, and developed my own code. See
http://www.mathworks.com/matlabcentral/newsreader/view_thread/296180
for details.
However, this may be not useful in this case since I also use factorization.

Subject: Alternatives/improvements to mldivide?

From: Svein Morten Drejer

Date: 7 Apr, 2011 07:25:19

Message: 10 of 10

"PLH" wrote in message <ibdn2c$j60$1@fred.mathworks.com>...
> Aha! In case it's useful for anyone else:
>
> http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6TYJ-4Y41V0R-3&_user=2195977&_coverDate=02%2F28%2F2010&_rdoc=1&_fmt=high&_orig=search&_origin=search&_sort=d&_docanchor=&view=c&_searchStrId=1534909065&_rerunOrigin=google&_acct=C000035218&_version=1&_urlVersion=0&_userid=2195977&md5=dd81b2a6736db70e1ace92cc0c850865&searchtype=a
>
> The above is a link to a paper detailing an algorithm designed to obtain the inverse (LU decomposition) for particular class of problem I stated above ("cyclic pentadiagonal" or "periodic pentadiagonal").
>
> PLH

Hi, I'm currently working on a very similar problem involving sparse matrices and the Crank-Nicholson method. The link to the paper appears to be broken, so I would be very gratefull if you could provide a working link.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us