mrdivide function in MATLAB: what is it doing, and how can I do it in Python?

Question

I have this line of MATLAB code:

a/b

I am using these inputs:

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

This is the result (a 1x25 matrix):

[5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]

What is MATLAB doing? I am trying to duplicate this behavior in Python, and the mrdivide documentation in MATLAB was unhelpful. Where does the 5 come from, and why are the rest of the values 0?

I have tried this with other inputs and receive similar results, usually just a different first element and zeros filling the remainder of the matrix. In Python when I use linalg.lstsq(b.T,a.T), all of the values in the first matrix returned (i.e. not the singular one) are 0.2. I have already tried right division in Python and it gives something completely off with the wrong dimensions.

I understand what a least square approximation is, I just need to know what mrdivide is doing.

Related:

• Array division- translating from MATLAB to Python

Answer 1

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

I did not find the earlier answers to create a satisfactory substitute, so I dug into Matlab's reference documents for mrdivide further and found the solution. I cannot explain the actual mathematics here or take credit for coming up with the answer. I'm just following Matlab's explanation. Additionally, I wanted to post the actual detail from Matlab to give credit. If it's a copyright issue, someone 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 that the ' symbol indicates the complex conjugate transpose. In python using numpy, that requires .conj().T chained together.

Answer 2

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

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

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

When the system is overdetermined, both algorithms provide the same answer. When the system is underdetermined, PINV will return the solution x, that has the minimum norm (min NORM(x)). MLDIVIDE will pick the solution with least number of non-zero 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 underdetermined, and the two different solutions 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 error of xb-a is non-negligible (non-exact solution) and the same, i.e. norm(x1*b-a) and norm(x2*b-a) will return the same result.

What is MATLAB doing?

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

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

What is Python doing?

Python, in linalg.lstsq uses pseudo-inverse/SVD, as demonstrated above (that's why you get a vector of 0.2's). In effect, the following will both give you 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))