I am trying to partially invert a non-square matrix M:V->W in the sense that, for a certain basis vector v \in V such that no other vector is mapped into Mv, I want to find a matrix N such that NMv = v. I say partial inversion because there may be other linearly independent vectors x, y such that Mx = My.
I have been using sympy to do this programmatically, and the only way I have found to do this is by using the pseudoinverse function .pinv(). However, this function hangs on a particular matrix that I want to (pseudo)invert -- I am unsure whether it is a bug or if the matrix is just too large.
Sympy can, however, reduce M to row echelon form with the .rref() function which runs very quickly. It would be nice to be able to extract the row operations (or elementary matrices) because they can easily be inverted to give the desired result. Is there some way to get the elementary matrices from the .rref() function (which also returns the pivots)? Is there another way unrelated to .pinv() to get the result I want?
I guess you mean something like this:
In [32]: x, y = symbols('x, y')
In [33]: M = Matrix([[1, 1], [0, 1], [0, 0]])
In [34]: xy = Matrix([x, y])
In [35]: XY = Matrix([-1, 1])
In [36]: M
Out[36]:
Matrix([
[1, 1],
[0, 1],
[0, 0]])
In [37]: M*xy
Out[37]:
Matrix([
[x + y],
[ y],
[ 0]])
In [38]: solve(Eq(M*xy, M*XY), [x, y])
Out[38]: {x: -1, y: 1}
You can find the matrix by forming the augmented matrix and using rref:
In [39]: Matrix.hstack(M, eye(3)).rref()
Out[39]:
⎛⎡1 0 1 -1 0⎤ ⎞
⎜⎢ ⎥ ⎟
⎜⎢0 1 0 1 0⎥, (0, 1, 4)⎟
⎜⎢ ⎥ ⎟
⎝⎣0 0 0 0 1⎦ ⎠
In [40]: Matrix.hstack(M, eye(3)).rref()[0][:2,2:]
Out[40]:
Matrix([
[1, -1, 0],
[0, 1, 0]])
In [41]: M.pinv()
Out[41]:
Matrix([
[1, -1, 0],
[0, 1, 0]])