Hi,

I am doing my master thesis at The University of Linköping (Sweden). Roughly: I am trying to create a direction field repressenting a mesh - and for that I need to solve a linear system Ax = b, i.e. I want to find x that solves the equation. Also the l2-norm of x must be minimized. The matrix A is sparse and underconstrained.

Is there anyway I can solve this in a elegant way, preferably without using the external solvers?

Best regards
Douglas Nordwall

If you known that the rows of A are linearly independent (so you have more columns than rows), then you can solve AA^T y = b (using any of our sparse linear solver) and get x = A^T y. Otherwise, you can use our SparseQR (in Eigen 3.2) which is rank revealing to help you computing the rank of A and assemble a full rank Gram matrix (AA^T here).

Are you sure this will work? When I try solving AA^Ty = b in matlab I get a vector of infinite valued numbers. I just want to make sure I get the math correct before trying to implement anything.

Code: Select all

A = [0.759836 0.759836 0 0.759836 0 0 ;
-0.759836 0 0.759836 0 0 0.759836 ;
0 -0.759836 -0.759836 0 0.759836 0 ;
0 0 0 -0.759836 -0.759836 -0.759836 ;];
 
b = [-3.14159;-3.14159;-3.14159;9.42478;];

y = A*A'\b

y =

   Inf
   Inf
   Inf
   Inf



Thank you for your time!
//Douglas

This is because rank(A)<rows(A), therefore A*A' is singular and you're back the original problem. How large are your matrices? If not too large, then a dense SVD will do the job easily. Otherwise, the best would be to perform a complete orthogonalization which is more tricky.

Hey, thanks for your answer! (Although my linear algebra is a bit outdated)

This is basicly what I wish to do:

[Q,R,E] = qr(A') ;
x = Q*(R'\(E'*b)) ;

Can I obtain Q, R and E with a SparseQR implementation?

//Douglas

Found Q, R, E using ColPivHouseholderQR.

Could you please elaborate the dense SVD solution? I am currently working on pretty small matrixes where M*N < 100. Im trying the JacobianSVD solver, and it doesnt give me the solution Im looking for.

The following matlab code produces the result I want:
[Q, R, E] = qr(A');
x = Q*(R'\(E'*b));

And I've managed to get the correct Q, R and E using Eigen.

Thanks for your time!
//Douglas

hm, which Eigen version are you using? JacobiSVD::solve should give you the minimal norm solution. Please try the latest Eigen version.

Regarding x = Q*(R'\(E'*b)); you can apply Q and E' using operator*. However, R'\ will require one more LS solver, e.g., using a QR factorization of R(1:rank,:)' which is full rank. A smarter approach is to eliminate R(1:rank,rank:columns) by applying householder reflectors on the right to get the full orthogonalization...

Yes, I have the latest release!

This is how I use the JacobiSVD:

Code: Select all

   
      //Setup A matrix   
      Eigen::MatrixXd Ae(AVector.size(), AVector[0].size());
      for(int i = 0; i < AVector.size(); i++)
      {
         for(int j = 0; j < AVector[0].size(); j++)
         {
            Ae(i, j) = AVector[i][j];
         }
      }

      //Setup b vector
      Eigen::VectorXd be(bVector.size()), xe;
      for(int i = 0; i < bVector.size(); i++)
         be(i) = bVector[i];

      //Solve x = A\b
      Eigen::JacobiSVD<Eigen::MatrixXd> svd(Ae, Eigen::ComputeFullU| Eigen::ComputeFullV);
      xe = svd.solve(be);

The vector xe contains a solution, but I dont think it is the minimum l2 norm solution. Am I using it wrong?


Ive also tried to QR factorize it first to solve x = Q*(R'\(E'*b)).

Code: Select all

      //x = Q*(R'\(E'*b)) ;
      Ae.transposeInPlace();
      Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(Ae);
      Q = qr.householderQ();
      R = qr.matrixQR().triangularView<Eigen::Upper>();
      E = qr.colsPermutation();

      E.transposeInPlace();
      R.transposeInPlace();
      
      E = E*be;
      Eigen::JacobiSVD<Eigen::MatrixXd> svd2(R, Eigen::ComputeThinU| Eigen::ComputeThinV);
      xe = Q*svd2.solve(E);
      cout << "xe: " << xe << "\n";


But without to much of a surprise it produces the same solution.

//Douglas

hm, our bad:

https://bitbucket.org/eigen/eigen/commits/99289da182ea/
Changeset: 99289da182ea
User: ggael
Date: 2013-11-01 18:21:46
Summary: Add a rank method with threshold control to JacobiSVD, and make solve uses it to return the minimal norm solution for rank-deficient problems

Works fine now! Thanks alot!