How can I reuse the SVD of matrix A to solve LS problems for both A and its transpose via Eigen C++?

If $A\in R^{m\times n}, b\in R^m, c\in R^n$,

if I need to solve the least square problems via SVD of $A$ and $A^T$, i.e.

I need to solve the least square solutions to following linear systems via SVD:

$Ax=b_i,\quad i=1,2,3,\cdots$

$A^Ty=c_j,\;j=1,2,3,\cdots$

it is also possible to do these via householderQr().solve() in Eigen, but it is instable; So I want to do this via SVD decomposition of $A$ and $A^T$; obviously, the SVD of $A^T$ can be obtained by the transpose of that for $A$.

My question is, how can I implement it in Eigen so that only SVD on $A$ is Ok to solve both the two type least square problems?

int m=100,n=50;
VectorXd b=VectorXd::randome(m),c=VectorXd::randome(n);
MatrixXd A=MatrixXd::random(m,n);
Eigen::JacobiSVD<Eigen::MatrixXd> _svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
MatrixXd  x = _svd.solve(b);

A.transposeInPlace();


// How can I save this sentence below which implements the SVD of $A^T$

Eigen::JacobiSVD<Eigen::MatrixXd> _svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);

MatrixXd  y = _svd.solve(c);

• This question is rather difficult to understand, since you didn't give any background. Are you asking how to solve a least squares problem using an SVD of A, or do you know how to do that and want help with doing that in Eigen? – Christian Clason Feb 6 '14 at 19:26

I do not have too much experience using eigen, but when you are solving the systems with the SVD decomposition, actually you are doing the followin:

$A x = USV^T x = b$

and you use the SVD decomposition of A to isolate x

$x = VS^{-1}U^T b$.

The same for the other system:

$A^T y = VSU^T y = c$

and you use the SVD decomposition of A to isolate x

$y = US^{-1}V^T c$.

So with eigen you should be able to recover a pointer to the matrices $U$, $S$ and $V$ as something similar to that

MatrixUType* U        = _svd.matrixU()
MatrixVType* V        = _svd.matrixV()
SingularValuesType* S = _svd.singularValues()


And then use some functions from the library to apply the operations to $c$ (matrix multiplications, traspose multimplication, and diagonal inversion) to obtain $y$.

I am sorry I can not give you more implementation details, but it think that it can help.

• thank you very much. I am just curious whether there is any solve() options something like applytoleft, applytoright so that it can solve the linear system of $A$, $A^T$. – LCFactorization Feb 8 '14 at 12:54