'Solve Ax = B using the Moore-Penrose pseudoinverse.
There may be zero, one, or infinite solutions.  If one solution
exists, it will be returned.  If infinite solutions exist, one will
be returned based on the value of arbitrary_matrix.  If no solutions
exist, the least-squares solution is returned.
Parameters
B : Matrix
The right hand side of the equation to be solved for.  Must have
the same number of rows as matrix A.
arbitrary_matrix : Matrix
If the system is underdetermined (e.g. A has more columns than
rows), infinite solutions are possible, in terms of an arbitrary
matrix.  This parameter may be set to a specific matrix to use
for that purpose; if so, it must be the same shape as x, with as
many rows as matrix A has columns, and as many columns as matrix
B.  If left as None, an appropriate matrix containing dummy
symbols in the form of ``wn_m`` will be used, with n and m being
row and column position of each symbol.
Returns
x : Matrix
The matrix that will satisfy Ax = B.  Will have as many rows as
matrix A has columns, and as many columns as matrix B.
Examples
>>> from sympy import Matrix
>>> A = Matrix([[1, 2, 3], [4, 5, 6]])
>>> B = Matrix([7, 8])
>>> A.pinv_solve(B)
Matrix([
[ _w0_0/6 - _w1_0/3 + _w2_0/6 - 55/18],
[-_w0_0/3 + 2*_w1_0/3 - _w2_0/3 + 1/9],
[ _w0_0/6 - _w1_0/3 + _w2_0/6 + 59/18]])
>>> A.pinv_solve(B, arbitrary_matrix=Matrix([0, 0, 0]))
Matrix([
[-55/18],
[   1/9],
[ 59/18]])
See Also
lower_triangular_solve
upper_triangular_solve
gauss_jordan_solve
cholesky_solve
diagonal_solve
LDLsolve
LUsolve
QRsolve
pinv
Notes
This may return either exact solutions or least squares solutions.
To determine which, check ``A * A.pinv() * B == B``.  It will be
True if exact solutions exist, and False if only a least-squares
solution exists.  Be aware that the left hand side of that equation
may need to be simplified to correctly compare to the right hand
side.
References
.. [1] https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#Obtaining_all_solutions_of_a_linear_system'

def pinv_solve(self, B, arbitrary_matrix=None):

    from sympy.matrices import eye
    A = self
    A_pinv = self.pinv()
    if (arbitrary_matrix is None):
        (rows, cols) = (A.cols, B.cols)
        w = symbols('w:{0}_:{1}'.format(rows, cols), cls=Dummy)
        arbitrary_matrix = self.__class__(cols, rows, w).T
    return ((A_pinv * B) + ((eye(A.cols) - (A_pinv * A)) * arbitrary_matrix))
