mrdivide function in Matlab- what is it doing and how to do it in Python?

Question

I have this line of 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: [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], a 1x25 matrix.

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 approx. is, I just need to know what mrdivide is doing.

Related:

• Array division- translating from MATLAB to Python

Answer 1

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

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

You can read more about least squares on wikipedia.

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

Answer 2

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

Answer 3

The rref() method of a sympy Matrix gets the Matlab result programmatically:

import sympy 

a = hstack([r_[1:10],r_[1:10]])
b = ones(25,18)
#use sympy's rref()
b2 = array(sympy.Matrix(b).rref()[0]).T
c = linalg.lstsq(b2,a)[0]

print c

[ 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.]

Since b is row rank 1, rref() keeps only one of the 25 rows. The rest become 0s. I transposed this to make it agree with the equation b2 * c = a, where 'b2' is size (18,25), 'c' is (25,) and a is (18,).

There are, of course, an infinite number of min norm solutions for this contrived problem. The mean output value is 5 (i.e. sum(a)/len(a)). Thus, given a matrix that's a row of ones, any vector whose components sum to 5 is a min norm solution. The two simplest answers are 5,0{24} and 0.2{25}.