Compute $\mathbf v \mathbf A^{-1}\mathbf v^\top$ in a numerically stable way

Question

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?

Answer 1

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).

• dot docs
• lstsq docs
• transpose docs

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

• cholesky

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.

Answer 2

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$$

Answer 3

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.)