13

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?
MattS
  • 138
  • 9
Schrodinger
  • 190
  • 1
  • 1
  • 7
  • 2
    What do you mean by "special solution"? The solution from python (`[0.918 3.541 1.279]`) is also a correct solution. You have 2 equations in 3 unkowns so there is no unique solution. The solutions are `[-1 9/2 0]+s*[3/2 -3/4 1]` for any real number `s`. Set `s=2/3` for the Matlab solution, and `s=1.27868852` for the Python solution. – David Nov 06 '15 at 05:23
  • Its not just about giving any correct solution. mldivide always gives a solution vector of n non-zero elements where n is the rank of the matrix, where as numpy doesn't do this. I'm looking to get the exact same solution in numpy as by mldivide. – Schrodinger Nov 06 '15 at 05:44
  • 2
    The Octave `A\b` solution is the same as the `numpy` one. MATLAB doc suggests `pinv(A)*B` as a computationally more expensive method. In Octave that produces the same thing. `numpy` also has `pinv`, witth the same result. – hpaulj Nov 06 '15 at 06:24
  • The operation you're looking at doesn't have a single numerical solution, so I don't see why you expect to obtain the same result if you ask various implementations to give a single answer to this more general question. (Although if you ask it right, MATLAB should also return a least-squares guess to the question) – Andras Deak -- Слава Україні Nov 06 '15 at 08:58
  • 4
    There is no least-squares solution to this---there are an infinite number of exact solutions. If something else in your code depends on getting the "right" solution, then you either need to specify another condition or there is a problem with your algorithm. – David Nov 06 '15 at 09:11
  • Oops, thanks @David, you're right. I confused it with the over-constrained case, mainly due to OP's use of `lstsq`...sorry. – Andras Deak -- Слава Україні Nov 06 '15 at 09:30
  • 2
    @David When there is an infinite number of exact solutions, `lstsq` returns the solution of minimal norm; i.e., one with the smallest sum of squares of the variables. This is something one might describe as "least squares", although "minimum-norm" is a preferable term to avoid confusion. –  Jul 06 '16 at 15:48
  • 1
    @Bookend I don't think I agree completely. A lest-squares solution means that the residual (error) in the approximate solution is minimised in a least-squares sense, i.e. the distance from the true solution to the approximate solution is minimised. Picking the solution with the smallest norm is a way to determine which solution to pick, just like Matlab tries to make as many elements of the vector zero as possible. So I think your method gives a least-squares solution, but uses a different way of picking which solution to return. – David Jul 06 '16 at 16:38

1 Answers1

15

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]]
  • 1
    It seems that `np.linalg.solve(mA, vB)` works only for squared matrices. Is there a solver for any size of matrix (Of course which has the same number of rows as `vB`)? – Royi Jun 25 '18 at 15:21