Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

I've read that you should avoid computing a matrix inverse, as you generally don't need to, but I don't know the best way to avoid it. I need to compute:

$$x = \mathbf v \mathbf A^{-1}\mathbf v^\top$$

where $x$ is a scalar, $\mathbf v$ is a row vector, $\mathbf A$ is a symmetric positive definite matrix (but perhaps with eigenvalues close to $0$) and ${}^\top$ means transpose.

I'm using numpy/scipy so feel free to express an answer using their functions.


EDIT:

Any pros/cons of the least squares approach versus doing an eigenvector decomposition?

share|improve this question
1  
The "safe" version of the least squares approach (SVD) and the eigendecomposition are the same, as I noted in percusse's answer. The "fast" version (QR, or Cholesky) is faster ;) , but a bit riskier if as you say your matrix has potentially tiny eigenvalues... –  J. M. Nov 23 '11 at 5:05
    
@J.M., thanks, speed is not an issue with me so I will go the safe route. –  Jeff Nov 23 '11 at 14:00

3 Answers 3

up vote 2 down vote accepted

v*inv(A) is the same as v/A in matlab notation, which uses a linear (least squares) solver rather than calculating the inverse.

So I guess I would use:

numpy.dot( numpy.linalg.lstsq( A, v ), v )

I'm not sure if the row versus column will matter to the software, but you can freely transpose things without changing anything important as A is symmetric (so if it cannot detect you passing the wrong kind of vector in, then use numpy.transpose).

If you have many different v for a single A, then you can use:

to get a better solver. I didn't see any such solver built-in to numpy, but things like matlab (and hopefully octave) would automatically switch to a faster solver for the SPD case.

share|improve this answer
    
I found cholesky in the parallel directory: decompositions. –  rcollyer Nov 23 '11 at 4:55

If the computational cost is not an issue, you can find the singular value decomposition of $A=U\Sigma V^T$ and only invert $\Sigma$ then do the multiplication as $$x=vU\Sigma^{-1}(Vv)^T$$

share|improve this answer
1  
The $\Sigma^{-1}$ can either be a true inverse or a pseudoinverse (reciprocate those diagonal entries larger than machine epsilon times the matrix norm), of course. Then again, since $A$ is symmetric positive definite, the SVD and the eigendecomposition are one and the same. –  J. M. Nov 23 '11 at 4:30

I've read that you should avoid computing a matrix inverse, as you generally don't need to.

This statement is incomplete. In general, you have no need for the inverse of a matrix as a separate entity. In most cases, you have something of the form $A^{-1}x$, where $x$ is a column vector, and computing $A^{-1}$ is as computationally intensive as computing $A^{-1}x$. So, if you don't need the inverse by itself, it is not worth calculating it in isolation.

Now when computing $A^{-1}x$, your actually computing the solution to $Ax = b$, and all numerical algorithms treat it as such. So, you need to use numpy.linalg.solve, as follows

x = numpy.linalg.solve( A, v )
numpy.dot( v, x )

Or, in one line

numpy.dot( v, numpy.linalg.solve( A, v ) )

This differs from Jack's solution, in that numpy.linalg.lstsq solves $A x = b$ by minimizing $|b - A x|^2_2$, while numpy.linalg.solve uses an LU decomposition and it assumes that $A$ has full rank. (Being symmetric implies that it is square.)

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.