python - Numpy vs mldivide,"" matlab operator

Question

A B in matlab gives a special solution while numpy.linalg.lstsq doesn't.

A = [1 2 0; 0 4 3];
b = [8; 18];
c_mldivide = A  b
c_mldivide =

                 0
                 4
  0.66666666666667
 c_lstsq = np.linalg.lstsq([[1 ,2, 0],[0, 4, 3]],[[8],[18]])
 print c_lstsq
 c_lstsq = (array([[ 0.91803279],
                   [ 3.54098361],
                   [ 1.27868852]]), array([], dtype=float64), 2, array([ 5.27316304,1.48113184]))

1. How does mldivide A B in matlab give a special solution?
2. Is this solution usefull in achieving computational accuracy?
3. Why is this solution special and how might you implement it in numpy?

Answer 1

For under-determined systems such as yours (rank is less than the number of variables), mldivide returns a solution with as many zero values as possible. Which of the variables will be set to zero is up to its arbitrary choice.

In contrast, the lstsq method returns the solution of minimal norm in such cases: that is, among the infinite family of exact solutions it will pick the one that has the smallest sum of squares of the variables.

So, the "special" solution of Matlab is somewhat arbitrary: one can set any of the three variables to zero in this problem. The solution given by NumPy is in fact more special: there is a unique minimal-norm solution

Which solution is better for your purpose depends on what your purpose is. The non-uniqueness of solution is usually a reason to rethink your approach to the equations. But since you asked, here is NumPy code that produces Matlab-type solutions.

import numpy as np
from itertools import combinations
A = np.matrix([[1 ,2, 0],[0, 4, 3]])
b = np.matrix([[8],[18]])

num_vars = A.shape[1]
rank = np.linalg.matrix_rank(A)
if rank == num_vars:              
    sol = np.linalg.lstsq(A, b)[0]    # not under-determined
else:
    for nz in combinations(range(num_vars), rank):    # the variables not set to zero
        try: 
            sol = np.zeros((num_vars, 1))  
            sol[nz, :] = np.asarray(np.linalg.solve(A[:, nz], b))
            print(sol)
        except np.linalg.LinAlgError:     
            pass                    # picked bad variables, can't solve

For your example it outputs three "special" solutions, the last of which is what Matlab chooses.

[[-1. ]
 [ 4.5]
 [ 0. ]]

[[ 8.]
 [ 0.]
 [ 6.]]

[[ 0.        ]
 [ 4.        ]
 [ 0.66666667]]