Thread Subject: large-scale sparse least square

Hi all,

I'm a student working in the field of computer vision.
I have a large -scale least square as follow:
x = M\q where M is nxn sparse banded , and n is 40,000 ~100,000.

Two questions:
1) Is "\" back slash the best method in matlab to solve this?
   or there is any iterative method can solve this
approximately which use less memory?

2) Forming M is painful and slow in matlab. But I have to
form M several times since the least square is a inner step
of Expectation Maximization algorithm. I might like to write
a Mex file to build M and solve x. Is there any source in
C++ about sparse matrices and large scale least square?

Best Regards,
Min Sun

"Min " <aliensunmin@gmail.com> wrote in message
<fjicau$6s$1@fred.mathworks.com>...

> or there is any iterative method can solve this
> approximately which use less memory?
>

Look at gmres, bicg, and ilu.

Bruno

I forgot to mention that if your problem is least-square
(linear model I suppose, since you invert matrix), it's
equivalent to solving a symmetric (auto-adjoint) problem.

A*x = y, lsq sense
<=> H*x = z, with H=A'*A, and z=A'y

H is symmetric definite positive.

You could use the well-known conjuguate gradient method to
solve the second formulation. Do not neglect preconditioning
in the large scale optimization.

Bruno

"Bruno Luong" <brunoluong@yahoo.com> wrote in message
<fjire8$9mp$1@fred.mathworks.com>...
> I forgot to mention that if your problem is least-square
> (linear model I suppose, since you invert matrix), it's
> equivalent to solving a symmetric (auto-adjoint) problem.
>
> A*x = y, lsq sense
> <=> H*x = z, with H=A'*A, and z=A'y
>
> H is symmetric definite positive.
>
> You could use the well-known conjuguate gradient method to
> solve the second formulation. Do not neglect preconditioning
> in the large scale optimization.
>
> Bruno

And I still forget one more remark. The H matrix should not
be explicitely computed (it could become much denser).
Because cgm only requires user to provide the calculation of
the product H*x, this can be accoplimshed by doing A'*(A*x).
This way uses best the sparseness of A.

Bruno

"Min " <aliensunmin@gmail.com> wrote in message
<fjicau$6s$1@fred.mathworks.com>...
> Hi all,
>
> I'm a student working in the field of computer vision.
> I have a large -scale least square as follow:
> x = M\q where M is nxn sparse banded , and n is 40,000
~100,000.
>
> Two questions:
> 1) Is "\" back slash the best method in matlab to solve this?
> or there is any iterative method can solve this
> approximately which use less memory?
>
> 2) Forming M is painful and slow in matlab. But I have to
> form M several times since the least square is a inner step
> of Expectation Maximization algorithm. I might like to write
> a Mex file to build M and solve x. Is there any source in
> C++ about sparse matrices and large scale least square?
>
> Best Regards,
> Min Sun


Forming M is slow only if you do it wrong. Are you using
statements such as

M(i,j)=x

? If so, then yes, it will be slow ... but that's the wrong
method. See Loren's March 1st blog on creating sparse matrices.

Creating M should be very fast unless you're doing it wrong.

On Dec 10, 12:55 pm, "Bruno Luong" <brunolu...@yahoo.com> wrote:
> I forgot to mention that if your problem is least-square
> (linear model I suppose, since you invert matrix), it's
> equivalent to solving a symmetric (auto-adjoint) problem.
>
> A*x = y, lsq sense
> <=> H*x = z, with H=A'*A, and z=A'y
>
> H is symmetric definite positive.
>
> You could use the well-known conjuguate gradient method to
> solve the second formulation. Do not neglect preconditioning
> in the large scale optimization.
>
> Bruno

YOU CAN USE NYSTROM APROXIMATION WITH SYMMLQ

You have TLSQR available in the TOMLAB Base Module:
http://tomopt.com/tomlab/products/base/solvers/

Check if spdiags can speed things up for you.

Best wishes, Marcus
Tomlab Optimization Inc.
http://tomopt.com/tomlab/

"Bruno Luong" <b.luong@fogale.fr> wrote in message
<fjj19c$21r$1@fred.mathworks.com>...
> "Bruno Luong" <brunoluong@yahoo.com> wrote in message
> <fjire8$9mp$1@fred.mathworks.com>...
> > I forgot to mention that if your problem is least-square
> > (linear model I suppose, since you invert matrix), it's
> > equivalent to solving a symmetric (auto-adjoint) problem.
> >
> > A*x = y, lsq sense
> > <=> H*x = z, with H=A'*A, and z=A'y
> >
> > H is symmetric definite positive.
> >
> > You could use the well-known conjuguate gradient method to
> > solve the second formulation. Do not neglect preconditioning
> > in the large scale optimization.
> >
> > Bruno
>
> And I still forget one more remark. The H matrix should not
> be explicitely computed (it could become much denser).
> Because cgm only requires user to provide the calculation of
> the product H*x, this can be accoplimshed by doing A'*(A*x).
> This way uses best the sparseness of A.
>

 Thanks a lot.
 My original problem is min x'Ax+2q'x
where A in theory should be Positive definite. However, when
the A is formed it's smallest eigenvalue is -1e-13.
What I previous did is solbe x = A\q anyway but since A has
very bad condition. x reture with lots of Nan and Inf.

I'm not very familiar with cg method. But I tried H = A'*A.
And H also have negative eigenvalue. I believe there are
some numerical error. Can cg handle this kind of problems?

Best
-MS

> Bruno

"Min " <aliensunmin@gmail.com> wrote in message
<fjkhlj$6f3$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.fr> wrote in message
> <fjj19c$21r$1@fred.mathworks.com>...
> > "Bruno Luong" <brunoluong@yahoo.com> wrote in message
> > <fjire8$9mp$1@fred.mathworks.com>...
> > > I forgot to mention that if your problem is least-square
> > > (linear model I suppose, since you invert matrix), it's
> > > equivalent to solving a symmetric (auto-adjoint) problem.
> > >
> > > A*x = y, lsq sense
> > > <=> H*x = z, with H=A'*A, and z=A'y
> > >
> > > H is symmetric definite positive.
> > >
> > > You could use the well-known conjuguate gradient method to
> > > solve the second formulation. Do not neglect
preconditioning
> > > in the large scale optimization.
> > >
> > > Bruno
> >
> > And I still forget one more remark. The H matrix should not
> > be explicitely computed (it could become much denser).
> > Because cgm only requires user to provide the calculation of
> > the product H*x, this can be accoplimshed by doing A'*(A*x).
> > This way uses best the sparseness of A.
> >
>
> Thanks a lot.
> My original problem is min x'Ax+2q'x
> where A in theory should be Positive definite. However, when
> the A is formed it's smallest eigenvalue is -1e-13.
> What I previous did is solbe x = A\q anyway but since A has
> very bad condition. x reture with lots of Nan and Inf.
>
> I'm not very familiar with cg method. But I tried H = A'*A.
> And H also have negative eigenvalue. I believe there are
> some numerical error. Can cg handle this kind of problems?
>
> Best
> -MS
>
> > Bruno
>
Thanks for all the helps.
I'm using pcg right now.
-MS

Can you send me a few of your matrices? I could use some
more computer-vision related problems in my collection:

http://www.cise.ufl.edu/research/sparse/matrices

You can upload the matrices to

http://www.cise.ufl.edu/~web-gfs/

Thanks,
Tim Davis (author of x=A\b when A is sparse and square)