[ACCEPTED]-mrdivide function in MATLAB: what is it doing, and how can I do it in Python?-linear-algebra

Accepted answer
Score: 13

MRDIVIDE or the / operator actually solves the xb = a linear 33 system, as opposed to MLDIVIDE or the \ operator 32 which will solve the system bx = a.

To solve a 31 system xb = a with a non-symmetric, non-invertible 30 matrix b, you can either rely on mridivide(), which 29 is done via factorization of b with Gauss 28 elimination, or pinv(), which is done via Singular 27 Value Decomposition, and zero-ing of the 26 singular values below a (default) tolerance 25 level.

Here is the difference (for the 24 case of mldivide): What is the difference between PINV and MLDIVIDE when I solve A*x=b?

When the system is overdetermined, both 23 algorithms provide the same answer. When 22 the system is underdetermined, PINV will 21 return the solution x, that has the minimum 20 norm (min NORM(x)). MLDIVIDE will pick 19 the solution with least number of non-zero 18 elements.

In your example:

% solve xb = a
a = [1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9];
b = ones(25, 18);

the system is 17 underdetermined, and the two different solutions 16 will be:

x1 = a/b; % MRDIVIDE: sparsest solution (min L0 norm) 
x2 = a*pinv(b); % PINV: minimum norm solution (min L2)

>> x1 = a/b
Warning: Rank deficient, rank = 1,  tol = 2.3551e-014.
ans =

    5.0000 0 0 ... 0 

>> x2 = a*pinv(b)
ans =

    0.2 0.2 0.2 ... 0.2 

In both cases the approximation 15 error of xb-a is non-negligible (non-exact solution) and 14 the same, i.e. norm(x1*b-a) and norm(x2*b-a) will return the same 13 result.

What is MATLAB doing?

A great break-down of the algorithms 12 (and checks on properties) invoked by the 11 '\' operator, depending upon the structure 10 of matrix b is given in this post in scicomp.stackexchange.com. I 9 am assuming similar options apply for the 8 / operator.

For your example, MATLAB is most 7 probably doing a Gaussian elimination, giving 6 the sparsest solution amongst a infinitude 5 (that's where the 5 comes from).

What is Python doing?

Python, in 4 linalg.lstsq uses pseudo-inverse/SVD, as demonstrated 3 above (that's why you get a vector of 0.2's). In 2 effect, the following will both give you 1 the same result as MATLAB's pinv():

from numpy import *

a = array([1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,8,9])
b = ones((25, 18))

# xb = a: solve b.T x.T = a.T instead 
x2 = linalg.lstsq(b.T, a.T)[0]
x2 = dot(a, linalg.pinv(b)) 
Score: 6

TL;DR: A/B = np.linalg.solve(B.conj().T, A.conj().T).conj().T

I did not find the earlier answers 12 to create a satisfactory substitute, so 11 I dug into Matlab's reference documents 10 for mrdivide further and found the solution. I 9 cannot explain the actual mathematics here 8 or take credit for coming up with the answer. I'm 7 just following Matlab's explanation. Additionally, I 6 wanted to post the actual detail from Matlab 5 to give credit. If it's a copyright issue, someone 4 tell me and I'll remove the actual text.

%/   Slash or right matrix divide.
%   A/B is the matrix division of B into A, which is roughly the
%   same as A*INV(B) , except it is computed in a different way.
%   More precisely, A/B = (B'\A')'. See MLDIVIDE for details.
%
%   C = MRDIVIDE(A,B) is called for the syntax 'A / B' when A or B is an
%   object.
%
%   See also MLDIVIDE, RDIVIDE, LDIVIDE.

%   Copyright 1984-2005 The MathWorks, Inc.

Note 3 that the ' symbol indicates the complex conjugate 2 transpose. In python using numpy, that requires 1 .conj().T chained together.

Score: 1

a/b finds the least square solution to the 11 system of linear equations bx = a

if b is 10 invertible, this is a*inv(b), but if it 9 isn't, the it is the x which minimises norm(bx-a)

You 8 can read more about least squares on wikipedia.

according 7 to matlab documentation, mrdivide will return at most k non-zero 6 values, where k is the computed rank of 5 b. my guess is that matlab in your case 4 solves the least squares problem given by 3 replacing b by b(:1) (which has the same 2 rank). In this case the moore-penrose inverse 1 b2 = b(1,:); inv(b2*b2')*b2*a' is defined and gives the same answer

Score: 1

Per this handy "cheat sheet" of numpy for matlab users, linalg.lstsq(b,a) -- linalg is 1 numpy.linalg.linalg, a light-weight version of the full scipy.linalg.

More Related questions