Using lsqnonlin on large-scale problems - Question about "JacobMult" option

Hi,

I am trying to use lsqnonlin() on a very large problem set, in which the Jacobian has a size of 31457280x393216 and is sparse. I am able to calculate the Jacobian myself and pass it to lsqnonlin() using

options = optimoptions('lsqnonlin','Jacobian','on')

The problem is that lsqnonlin() is taking a very long time. I stepped through the code to find where the bottleneck is happening and have tracked it down to the preconditioner step in the trust region algorithm. That is, within the function trdog(), it calls a preconditioner function aprecon() which performs a QR decomposition on the Jacobian. Specificially it is this section of aprecon() which takes an incredibly long time:

   p = colamd(TM);
   RPCMTX = qr(TM(:,p));

My question is, will using the JacobMult option of lsqnonlin() avoid this? According to the documentation:

 For large-scale structured problems, this function computes the Jacobian matrix product J*Y,   J'*Y, or J'*(J*Y) without actually forming J.

which to me sounds like it is only useful for large-scale problems in which the Jacobian is too large to even be stored, which is not the issue in my case. Additionally, I have read this: Jacobian Multiply Function with Linear Least Squares, and am confused on what the Jacobian Multiply function is actually doing, and therefore am unsure on how to even implement this.

To summarize, will the JacobMult option avoid the extremely time-consuming QR decomposition in lsqnonlin() for large-scale problems, and can anyone give a clear demonstration on what the JacobMult option actually is doing and how to implement it?

Thanks.

Answer 1

[Matt J]

If you use JacobMult, I imagine it will try to solve for the Newton steps iteratively, instead of with QR decomposition, which could be faster, depending on how sparse your Jacobian really is.

Basically, the Newton step amounts to solving a linear equation

J*x=F

where J is your Jacobian. There are iterative algorithms which can do this using a series of multiplications with J and J.'. In your case, the JacobMult function could be implemented simply as

function W=jmfun(J,Y,flag)

   if flag==0

    W = J'*(J*Y);

   elseif flag>0

     W= J*Y;

   else%flag<0

     W=J.'*Y;

   end

end

[Shayan]

Thanks, that clarifies things a bit. However after implementing this and looking through the documentation, it appears that lsqnonlin() always uses the first argument of the Jacobian-Multiply function to calculate a preconditioner using QR. Therefore if I simply use the Jacobian I originally calculated I run into the same problem. I think what I need to do is implement something like Jinfo without actually calculating the Jacobian. According to the documentation:

 Function handle for Jacobian multiply function. For large-scale structured problems, this function computes the Jacobian matrix product J*Y, J'*Y, or J'*(J*Y) without actually forming J. The function is of the form

 W = jmfun(Jinfo,Y,flag) 

 where Jinfo contains the matrix used to compute J*Y (or J'*Y, or J'*(J*Y))

 ...

 In each case, J is not formed explicitly. lsqnonlin uses Jinfo to compute the preconditioner.

Looking at the aprecon() code, if Jinfo is not equal to the Jacobian, but is instead something else with a different size (but still correctly calculates W in the Jacobian Multiply function), it will skip QR decomposition. But how do I choose a Jinfo that does this?

[Matt J]

That sounds weird, if it's true. I suppose you could trick aprecon by computing Jinfo as

   [m,n]=size(Jacobian);
   [i,j,s]=find(Jacobian);
   Jinfo={i,j,s,m,n};

and then rebuilding the Jacobian inside jmfun

   J=sparse(Jinfo{:});

   W=...

You could also perhaps just reshape() the Jacobian into a different sized matrix Jinfo. Undoing that in jmfun might be faster.

In any case, bear in mind that not using a pre-conditioner could end up slowing convergence to a solution, with no clear net gains in speed.

[Shayan]

aprecon() still calculates a preconditioner, but it appears to calculate a diagonal preconditioner based on the results of the Jacobian Multiply function instead of QR decomposition. Passing Jinfo as a cell produces an error since lsqnonlin() expects it to be of type double, but I tried out your idea of reshaping the Jacobian and undoing it in jmfun. It seems to work, the iterations are proceeding without getting stuck in QR decomposition.

Thanks!

Answer 2

[françois anquez]

Dear Shayan, Dear Matt,

I do not wether this will be usefull given the time past since you posted.

I was facing the same problem.

Something that helped for me was to set : 'SubproblemAlgorithm' to 'cg'.

This seems to considerably reduce the time spent in aprecon with similar results.

François.