I did a small experiment the other week and managed to invert some small matrices by the following procedure:
- Left hand side = matrix we want to invert (say $A$). If non-symmetric, instead $A^TA$
- First standard basis vector $v_1 = [1,0,\cdots,0]^T$ becomes RHS if non-symmetric LHS, if symmetric then $A^Tv_1$ instead.
- Solve with conjugate-gradient.
- Rotate the $1$ in rhs one spot $v_2=[0,1,0,\cdots,0]^T$ becomes RHS if non-symmetric LHS, if symmetric then $A^Tv_2$
- Repeat for all columns of standard basis vectors: Solve with Conjugate-Gradient.
- Collect the solution vectors as columns in the new matrix.
- Done!
Everything will be extremely parallellizable. Each column could be assigned a set of processing units because each column is independent of the solvin of the other columns. Everything to write and accumulate will be local in memory for the (groups of) cores. The matrix to invert can be stored in common "global memory".
(Many modern hardwares, for example GPU architectures has this setup)
Do there already exist software packages which do this? If not, why not? Is there some obvious drawback that I am not seeing?
EDIT:
Example where it is interesting to solve for a subset of rows or columns of the inverse: We discretize a differential equation. The solution to the linear equation system for this discretized differential equation is the same as to multiply the inverse matrix with a vector (right-hand-side). Then one row in the inverse is the linear combination of "initial/boundary" conditions for the solution at that particular point in space and/or time. For different rows we can see how this varies throughout the definition space in a very illustrative way.
Here is a motivational gif for why it is useful. The equation system it illustrates the solution for is this Nyquist sampling question. We can see the kernel shifting to become the linear interpolation weight as the output position varies over the signal.
