Thread Subject: nonnegative Ax=b lsq for large, sparse A.

Hello all,

I have a large, sparse, rectangular (overdetermined) matrix A.

I'm solving Ax = b, where all elements of A, x and b are nonnegative.

I thought at first that lsqnonneg would be my ideal function, however, I hadn't realised that it is not useful for large sparse matrices (creates a full matrix the same size as A).

Does anyone know of an alternative for lsqnonneg which may be used with a large, sparse coefficients matrix?

A sample of the problem (.mat file containing A and b) can be found in:
http://www2.eng.cam.ac.uk/~thc29/lsqnonneg_sparse/

Tim Davis, if you happen across this, you may enjoy the sparsity pattern (it's for a tomographic reconstruction of a particle field in fluid flow).

Thanks all, and kind regards

Tom Clark

Incidentally, A x and b are real, too :)

- Tom C

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22iil$a5s$1@fred.mathworks.com>...
> Hello all,
>
> I have a large, sparse, rectangular (overdetermined) matrix A.
>
> I'm solving Ax = b, where all elements of A, x and b are nonnegative.
>
> I thought at first that lsqnonneg would be my ideal function, however, I hadn't realised that it is not useful for large sparse matrices (creates a full matrix the same size as A).
>
> Does anyone know of an alternative for lsqnonneg which may be used with a large, sparse coefficients matrix?
>
> A sample of the problem (.mat file containing A and b) can be found in:
> http://www2.eng.cam.ac.uk/~thc29/lsqnonneg_sparse/
>
> Tim Davis, if you happen across this, you may enjoy the sparsity pattern (it's for a tomographic reconstruction of a particle field in fluid flow).
>
> Thanks all, and kind regards
>
> Tom Clark

lsqlin does accept sparse bound constrained problems.

A = sprand(2000,1000,.005);
b = rand(2000,1);
lb = zeros(1000,1);

tic,x = A\b;toc
Elapsed time is 27.267806 seconds.

tic,x = lsqlin(A,b,[],[],[],[],lb,[],[],opts);toc
Optimization terminated: relative function value changing by less
 than sqrt(OPTIONS.TolFun), and rate of progress is slow.
Elapsed time is 2.186524 seconds.

John

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22iil$a5s$1@fred.mathworks.com>...
> Hello all,
>
> I have a large, sparse, rectangular (overdetermined) matrix A.
>
> I'm solving Ax = b, where all elements of A, x and b are nonnegative.
>
> I thought at first that lsqnonneg would be my ideal function, however, I hadn't realised that it is not useful for large sparse matrices (creates a full matrix the same size as A).
>
> Does anyone know of an alternative for lsqnonneg which may be used with a large, sparse coefficients matrix?
>
> A sample of the problem (.mat file containing A and b) can be found in:
> http://www2.eng.cam.ac.uk/~thc29/lsqnonneg_sparse/
>
> Tim Davis, if you happen across this, you may enjoy the sparsity pattern (it's for a tomographic reconstruction of a particle field in fluid flow).
>
> Thanks all, and kind regards
>
> Tom Clark

Gee, how did you know I'd be interested in large sparse matrices ... ? :-)

I did download it, but it didn't load properly. I got very bizarre error messages:
HDF5-DIAG: Error detected in HDF5 (1.8.1) thread 0:
  #000: H5Gdeprec.c line 1063 in H5Gget_objname_by_idx(): not a location ID
    major: Invalid arguments to routine
    minor: Inappropriate type
  #001: H5Gloc.c line 241 in H5G_loc(): invalid object ID
    major: Invalid arguments to routine
    minor: Bad value

And I get a blizzard of those errors (living in Florida, I guess I should say "hurricane" not "blizzard", but I grew up in Indiana...).

Can you do this instead?

[i j x] = find (A) ;
f=fopen('A.txt','w') ;
fprintf (f,'%d %d %32.17g\n', [i j x]') ;
fclose (f) ;
f=fopen('b.txt','w') ;
fprintf (f,'%32.17g\n',b) ;
fclose (f)

assuming b is a dense column vector, and post these txt files instead?

And may I add it to my collection?

John ... which version MATLAB are you using? I tried your example on a Pentium 4 running MATLAB 7.8 (R2009a) and sparse backslash took just 4.2 seconds. lsqlin took 0.9 seconds (with no "opts" argument since you didn't say what it was). SPQR_SOLVE took 0.6 seconds, but of course A\b and spqr_solve return x with negative values, whereas lsqlin respects the lower bound lb.

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22iil$a5s$1@fred.mathworks.com>...
> Hello all,
>
> I have a large, sparse, rectangular (overdetermined) matrix A.
>
> I'm solving Ax = b, where all elements of A, x and b are nonnegative.


Here's a very simple Majorize-Minimize approach. How fast it will be depends on the structure of A, of course

 Curvs=A'*(sum(A,2)); %majorizing curvatures

 x=x0; %Initialize x with some vector x0>=0

for i=1:NumIterations %MAIN LOOP

 gradient=A.'*(Ax-b);
 x=x-gradient./Curvs;
 x(x<0)=0;

end

"Tim Davis" <davis@cise.ufl.edu> wrote in message <h22o30$gi5$1@fred.mathworks.com>...
>
> John ... which version MATLAB are you using? I tried your example on a Pentium 4 running MATLAB 7.8 (R2009a) and sparse backslash took just 4.2 seconds. lsqlin took 0.9 seconds (with no "opts" argument since you didn't say what it was). SPQR_SOLVE took 0.6 seconds, but of course A\b and spqr_solve return x with negative values, whereas lsqlin respects the lower bound lb.

I'm usually running an older Mac, so R2007b is the
latest I can use, except when I can pry the laptop away
from the household CEO.

Worse, I was kibbitzing a bridge event on the bbo
vugraph at the time, so that was probably stealing a
few processor cycles away.

John

>"John D'Errico" <woodchips@rochester.rr.com> wrote in message
>*snip*
> lsqlin does accept sparse bound constrained problems.
>
> A = sprand(2000,1000,.005);
> b = rand(2000,1);
> lb = zeros(1000,1);
>
> tic,x = A\b;toc
> Elapsed time is 27.267806 seconds.
>
> tic,x = lsqlin(A,b,[],[],[],[],lb,[],[],opts);toc
> Optimization terminated: relative function value changing by less
> than sqrt(OPTIONS.TolFun), and rate of progress is slow.
> Elapsed time is 2.186524 seconds.
>
> John

John, thanks for this tip.

Unfortunately, lsqlin isn't quite there... The solution is pathologically slow for this matrix - my feeling is that the preconditioning isn't working well. I've tried increasing the PrecondBandWidth option, but to little avail.(increasing it to Inf causes a direct factorisation which leaves me out of memory, of course).

Has anyone got any more ideas?

Thanks

Tom Clark

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22ot5$bca$1@fred.mathworks.com>...

> Unfortunately, lsqlin isn't quite there... The solution is pathologically slow for this matrix - my feeling is that the preconditioning isn't working well. I've tried increasing the PrecondBandWidth option, but to little avail.(increasing it to Inf causes a direct factorisation which leaves me out of memory, of course).
--------

If the Hessian A.'*A is diagonally concentrated, you might be able to get away with pre-conditioning by inverting the Hessian diagonal. Here's a modification of the Maj-Min approach incorporating this idea.

W=sum(A.^2,2); %Hessian diagonal
W=1./W; %invert



 Curvs=A'*(sum(A,2).*W); %majorizing curvatures

 x=x0; %Initialize x with some vector x0>=0

for i=1:NumIterations %MAIN LOOP

 gradient=A.'*((Ax-b).*W);
 x=x-gradient./Curvs;
 x(x<0)=0;

end

@Tim,

Thanks very much, you're more than welcome to add it to your collection - if you like, I can submit it via the online interface.
I've resubmitted the matrix in text file form (A and b) at the same link, the .mat file must have become corrupted.
www2.eng.cam.ac.uk/~thc29/lsqnonneg_sparse/

@Matt,

Unfortunately, the hessian is not diagonally concentrated, so I'm trying the solution out without a preconditioner.

However, I am getting some very good results. This is standing up extremely well against an alternative (iterative) approach to the same overall problem... It's also reasonably quick even without preconditioning. I'm on the way, and will invest more time testing this rigorously.

Thanks, both of you :)

Tom Clark



"Matt " <xys@whatever.com> wrote in message <h22q2l$sdn$1@fred.mathworks.com>...
> "Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22ot5$bca$1@fred.mathworks.com>...
>
> > Unfortunately, lsqlin isn't quite there... The solution is pathologically slow for this matrix - my feeling is that the preconditioning isn't working well. I've tried increasing the PrecondBandWidth option, but to little avail.(increasing it to Inf causes a direct factorisation which leaves me out of memory, of course).
> --------
>
> If the Hessian A.'*A is diagonally concentrated, you might be able to get away with pre-conditioning by inverting the Hessian diagonal. Here's a modification of the Maj-Min approach incorporating this idea.
>
> W=sum(A.^2,2); %Hessian diagonal
> W=1./W; %invert
>
>
>
> Curvs=A'*(sum(A,2).*W); %majorizing curvatures
>
> x=x0; %Initialize x with some vector x0>=0
>
> for i=1:NumIterations %MAIN LOOP
>
> gradient=A.'*((Ax-b).*W);
> x=x-gradient./Curvs;
> x(x<0)=0;
>
> end

Tom,

lsqnonneg was upgraded in version R2008b; it now uses sparse linear algebra if the
matrix is sparse (it uses "\", backslash).

R2008b Release notes:
http://www.mathworks.com/access/helpdesk/help/techdoc/rn/brqyzsl-1.html
(See "Upgrades to lsqnonneg".)

-Marcelo

Thomas Clark wrote:
> Hello all,
>
> I have a large, sparse, rectangular (overdetermined) matrix A.
>
> I'm solving Ax = b, where all elements of A, x and b are nonnegative.
>
> I thought at first that lsqnonneg would be my ideal function, however, I hadn't realised that it is not useful for large sparse matrices (creates a full matrix the same size as A).
>
> Does anyone know of an alternative for lsqnonneg which may be used with a large, sparse coefficients matrix?
>
> A sample of the problem (.mat file containing A and b) can be found in:
> http://www2.eng.cam.ac.uk/~thc29/lsqnonneg_sparse/
>
> Tim Davis, if you happen across this, you may enjoy the sparsity pattern (it's for a tomographic reconstruction of a particle field in fluid flow).
>
> Thanks all, and kind regards
>
> Tom Clark

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22tdi$em2$1@fred.mathworks.com>...

> @Matt,
>
> Unfortunately, the hessian is not diagonally concentrated, so I'm trying the solution out without a preconditioner.
>
> However, I am getting some very good results. This is standing up extremely well against an alternative (iterative) approach to the same overall problem... It's also reasonably quick even without preconditioning. I'm on the way, and will invest more time testing this rigorously.
-----

That's good news, I guess, but also puzzling. The theory says that if the quadratic is both ill-conditioned and non-diagonally concentrated, the algorithm should be quite slow. Oh well...

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22tdi$em2$1@fred.mathworks.com>...
> @Tim,
>
> Thanks very much, you're more than welcome to add it to your collection - if you like, I can submit it via the online interface.
> I've resubmitted the matrix in text file form (A and b) at the same link, the .mat file must have become corrupted.
> www2.eng.cam.ac.uk/~thc29/lsqnonneg_sparse/
>
> @Matt,
>
> Unfortunately, the hessian is not diagonally concentrated, so I'm trying the solution out without a preconditioner.
>
> However, I am getting some very good results. This is standing up extremely well against an alternative (iterative) approach to the same overall problem... It's also reasonably quick even without preconditioning. I'm on the way, and will invest more time testing this rigorously.
>
> Thanks, both of you :)
>
> Tom Clark


I've downloaded the matrix (A.txt and b.txt). Both x=A\b and x=spqr_solve(A,b) ran out of memory on my 32-bit, 4GB RAM, 3.2GHz Pentium 4 computer with MATLAB R2009a. And of course they woudn't give the right solution anyway even if they did work. I then tried lsqlin (A,b, [],[],[],[], lb) ; where lb is an all-zeros vector. It finished in 88 seconds, and returns an x with norm(A*x-b) of about 29.97, and min(x) is about 1e-15.

I tried lsqnonneg(A,b), after waiting maybe 10 minutes, it's is still running. It would go faster if the backslash inside is replaced with spqr_solve. Maybe I'll give that a try.

Once you figure out a solution using Matt's method, could you post, here, a function

x = mysolver(A,b) ;

that solves your system? And maybe a flag saying whether or not to do preconditioning? Then I can post that along with the matrix, in my collection, and then if someone looks at the problem there they might find yet a better method. In any case, it's nice to include a proper solution method with a matrix in my collection.

Tim,

Thanks for your investigation. If you hang on as long as I did, you'll get frustrated with lsqnonneg.

My strategy will be to write a version of lsqnonneg ( I have the textbook containing the original 'nnls' algorithm on which it is based ), with the ability to use various different solution techniques from both the MATLAB tools and your toolbox(es). I think it'll make a nice addition to the FEX, and I'll trial the different methods on the matrix in question. Once I've done that, I'll let you know a solution method and preconditioning options.

I'll be away for a couple of weeks, so will post once I'm back.

Thanks, and kind regards

Tom Clark



 

"Tim Davis" <davis@cise.ufl.edu> wrote in message <h2go35$kl1$1@fred.mathworks.com>...
> "Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h22tdi$em2$1@fred.mathworks.com>...
> > @Tim,
> >
> > Thanks very much, you're more than welcome to add it to your collection - if you like, I can submit it via the online interface.
> > I've resubmitted the matrix in text file form (A and b) at the same link, the .mat file must have become corrupted.
> > www2.eng.cam.ac.uk/~thc29/lsqnonneg_sparse/
> >
> > @Matt,
> >
> > Unfortunately, the hessian is not diagonally concentrated, so I'm trying the solution out without a preconditioner.
> >
> > However, I am getting some very good results. This is standing up extremely well against an alternative (iterative) approach to the same overall problem... It's also reasonably quick even without preconditioning. I'm on the way, and will invest more time testing this rigorously.
> >
> > Thanks, both of you :)
> >
> > Tom Clark
>
>
> I've downloaded the matrix (A.txt and b.txt). Both x=A\b and x=spqr_solve(A,b) ran out of memory on my 32-bit, 4GB RAM, 3.2GHz Pentium 4 computer with MATLAB R2009a. And of course they woudn't give the right solution anyway even if they did work. I then tried lsqlin (A,b, [],[],[],[], lb) ; where lb is an all-zeros vector. It finished in 88 seconds, and returns an x with norm(A*x-b) of about 29.97, and min(x) is about 1e-15.
>
> I tried lsqnonneg(A,b), after waiting maybe 10 minutes, it's is still running. It would go faster if the backslash inside is replaced with spqr_solve. Maybe I'll give that a try.
>
> Once you figure out a solution using Matt's method, could you post, here, a function
>
> x = mysolver(A,b) ;
>
> that solves your system? And maybe a flag saying whether or not to do preconditioning? Then I can post that along with the matrix, in my collection, and then if someone looks at the problem there they might find yet a better method. In any case, it's nice to include a proper solution method with a matrix in my collection.

In article <h231jk$jcv$1@fred.mathworks.com>,
mREMOVEmaALLraCAPITALSzzi@mathworks.com says...
> Tom,
>
> lsqnonneg was upgraded in version R2008b; it now uses sparse linear algebra if the
> matrix is sparse (it uses "\", backslash).
>
> R2008b Release notes:
> http://www.mathworks.com/access/helpdesk/help/techdoc/rn/brqyzsl-1.html
> (See "Upgrades to lsqnonneg".)
>
> -Marcelo
>
Tom-

Can you tell us what version you are using? If earlier than 8b, can you
try it in 8b or later?


--
Loren
http://blogs.mathworks.com/loren

Loren,

It's 2008a. I'm trying to get my sysadmin to upgrade me to 9a so that I get the more recent version of lsqnonneg, amongst other things.

However, I have a serious feeling that even the newer version (which uses backslash so handles large sparse systems) will be too slow for me, because of the matrix structure (I need to repeat this tens of thousands of times) - hence I'm thinking that a more versatile version of nnls, which can call different preconditioners / solvers etc would be a useful tool - even shaving a few seconds off a solution is valuable.

Thanks!

Tom Clark


Loren Shure <loren@mathworks.com> wrote in message <MPG.24b661b2898ecd509899fb@news.mathworks.com>...
> In article <h231jk$jcv$1@fred.mathworks.com>,
> mREMOVEmaALLraCAPITALSzzi@mathworks.com says...
> > Tom,
> >
> > lsqnonneg was upgraded in version R2008b; it now uses sparse linear algebra if the
> > matrix is sparse (it uses "\", backslash).
> >
> > R2008b Release notes:
> > http://www.mathworks.com/access/helpdesk/help/techdoc/rn/brqyzsl-1.html
> > (See "Upgrades to lsqnonneg".)
> >
> > -Marcelo
> >
> Tom-
>
> Can you tell us what version you are using? If earlier than 8b, can you
> try it in 8b or later?
>
>
> --
> Loren
> http://blogs.mathworks.com/loren

Tom,

I have written a non-negative projected conjugate gradient solver that could be a suitable method fore sparse matrix. Unfortunately for various reasons I cannot release to source code. If you are interested I can investigate how well it performs. Not that the only thing I could provide is a p-code. If such solution is not acceptable, then I will not make any investigation.

Bruno

Bruno,

That's really generous - I think that pcg might be the way to go, so am planning on incorporating that.

Unfortunately, I'm building it into a larger toolbox (for tomographic particle image velocimetry) which I'll ultimately be distributing in academic circles, so I need source code. Also, I'm developing a new technique for flow measurement using this solution, and I'll need to explain my algorithm to the scientific community (and my PhD examiners!).

Probably best, then, to save your time rather than investigating your routine's properties on this matrix.

Thankyou very much, however - that's very kind.

Tom Clark




"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h2i7ac$89v$1@fred.mathworks.com>...
> Tom,
>
> I have written a non-negative projected conjugate gradient solver that could be a suitable method fore sparse matrix. Unfortunately for various reasons I cannot release to source code. If you are interested I can investigate how well it performs. Not that the only thing I could provide is a p-code. If such solution is not acceptable, then I will not make any investigation.
>
> Bruno

Well, this may not help you either, but I'll try anyways...

If you have the Optimization Toolbox, there is a linear least squares solver (lsqlin) that solves bound constrained problems using PCG. Granted, this solver is not dedicated to non-negative least-squares a la lsqnonneg or nnls, but it's likely to be able to handle larger problems.

-Steve

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h2ibjb$ht7$1@fred.mathworks.com>...
> Bruno,
>
> That's really generous - I think that pcg might be the way to go, so am planning on incorporating that.
>
> Unfortunately, I'm building it into a larger toolbox (for tomographic particle image velocimetry) which I'll ultimately be distributing in academic circles, so I need source code. Also, I'm developing a new technique for flow measurement using this solution, and I'll need to explain my algorithm to the scientific community (and my PhD examiners!).
>
> Probably best, then, to save your time rather than investigating your routine's properties on this matrix.
>
> Thankyou very much, however - that's very kind.
>
> Tom Clark
>
>
>
>
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h2i7ac$89v$1@fred.mathworks.com>...
> > Tom,
> >
> > I have written a non-negative projected conjugate gradient solver that could be a suitable method fore sparse matrix. Unfortunately for various reasons I cannot release to source code. If you are interested I can investigate how well it performs. Not that the only thing I could provide is a p-code. If such solution is not acceptable, then I will not make any investigation.
> >
> > Bruno

"Thomas Clark" <t.clark@remove.spamcantab.net> wrote in message <h2htl1$38n$1@fred.mathworks.com>...
> Tim,
>
> Thanks for your investigation. If you hang on as long as I did, you'll get frustrated with lsqnonneg.
>
> My strategy will be to write a version of lsqnonneg ( I have the textbook containing the original 'nnls' algorithm on which it is based ), with the ability to use various different solution techniques from both the MATLAB tools and your toolbox(es). I think it'll make a nice addition to the FEX, and I'll trial the different methods on the matrix in question. Once I've done that, I'll let you know a solution method and preconditioning options.
>
> I'll be away for a couple of weeks, so will post once I'm back.
>
> Thanks, and kind regards
>
> Tom Clark


lsqlin worked great for this problem.

lsqnonneg ran all night, in MATLAB R2009a. I'll try replacing backslash (which uses a slow MATLAB sparse QR) with x=spqr_solve(A,b), which is faster.

Failing that, I'll try my even spiffier MouseHolder QR:

http://www.cise.ufl.edu/~davis/Poetry/Mouseholder_QR.html

;-)