Array division- translating from MATLAB to Python

Go To StackoverFlow.com

7

I have this line of code in MATLAB, written by someone else:

c=a.'/b

I need to translate it into Python. a, b, and c are all arrays. The dimensions that I am currently using to test the code are:

a: 18x1,
b: 25x18,

which gives me c with dimensions 1x25.

The arrays are not square, but I would not want the code to fail if they were. Can someone explain exactly what this line is doing (mathematically), and how to do it in Python? (i.e., the equivalent for the built-in mrdivide function in MATLAB if it exists in Python?)

2009-06-16 13:54
by EmilyS
I think you have a typo. If "a" is 1-by-18, you don't need the transpose - gnovice 2009-06-16 14:13
It's not a typo, the Matlab code works perfectly - EmilyS 2009-06-16 14:31
@Emily: Then "a" has to be 18-by-1 (before the transpose), not 1-by-18. Otherwise MATLAB throws an error - gnovice 2009-06-16 14:36
Yes, A-transpose/B to rightly occur, A-transponse and B should have the same no. of columns. And so, A's size should be 18x1 as gnovice pointed out - Suvesh Pratapa 2009-06-16 14:43
Ok, sorry about that, it is 18x1 at first - EmilyS 2009-06-16 14:49


7

The line

c = a.' / b

computes the solution of the equation c b = aT for c. Numpy does not have an operator that does this directly. Instead you should solve bT cT = a for cT and transpose the result:

c = numpy.linalg.lstsq(b.T, a.T)[0].T
2009-06-17 18:41
by Vebjorn Ljosa
You may have the "a" and "b" in your equation flipped relative to the OPs "a" and "b". The line "c = a.'/b" is solving an equation of the form "xb = a.'", which becomes "(b.')(x.') = a". The OPs "b" (i.e. the matrix of equation coefficients) should be the first input to lstsq, transposed of course - gnovice 2009-06-17 19:20
Thanks for pointing that out. I'll edit my answer - Vebjorn Ljosa 2009-06-18 18:29


7

The symbol / is the matrix right division operator in MATLAB, which calls the mrdivide function. From the documentation, matrix right division is related to matrix left division in the following way:

B/A = (A'\B')'

If A is a square matrix, B/A is roughly equal to B*inv(A) (although it's computed in a different, more robust way). Otherwise, x = B/A is the solution in the least squares sense to the under- or over-determined system of equations x*A = B. More detail about the algorithms used for solving the system of equations is given here. Typically packages like LAPACK or BLAS are used under the hood.

The NumPy package for Python contains a routine lstsq for computing the least-squares solution to a system of equations. This routine will likely give you comparable results to using the mrdivide function in MATLAB, but it is unlikely to be exact. Any differences in the underlying algorithms used by each function will likely result in answers that differ slightly from one another (i.e. one may return a value of 1.0, whereas the other may return a value of 0.999). The relative size of this error could end up being larger, depending heavily on the specific system of equations you are solving.

To use lstsq, you may have to adjust your problem slightly. It appears that you want to solve an equation of the form cB = a, where B is 25-by-18, a is 1-by-18, and c is 1-by-25. Applying a transpose to both sides gives you the equation BTcT = aT, which is a more standard form (i.e. Ax = b). The arguments to lstsq should be (in this order) BT (an 18-by-25 array) and aT (an 18-element array). lstsq should return a 25-element array (cT).

Note: while NumPy doesn't make any distinction between a 1-by-N or N-by-1 array, MATLAB certainly does, and will yell at you if you don't use the proper one.

2009-06-16 14:07
by gnovice
Numpy does not make any distinction between a 1-by-N or N-by-1 array. A 1D array is a 1D array. If you make a 2D array with one dimension of length 1, then you can differentiate, but there's usually no reason to do so - endolith 2013-12-20 02:36


5

In Matlab, A.' means transposing the A matrix. So mathematically, what is achieved in the code is AT/B.


How to go about implementing matrix division in Python (or any language) (Note: Let's go over a simple division of the form A/B; for your example you would need to do AT first and then AT/B next, and it's pretty easy to do the transpose operation in Python |left-as-an-exercise :)|)

You have a matrix equation C*B=A (You want to find C as A/B)

RIGHT DIVISION (/) is as follows:

C*(B*BT)=A*BT

You then isolate C by inverting (B*BT)

i.e.,

C = A*BT*(B*BT)' ----- [1]

Therefore, to implement matrix division in Python (or any language), get the following three methods.

  • Matrix multiplication
  • Matrix transpose
  • Matrix inverse

Then apply them iteratively to achieve division as in [1].

Only, you need to do AT/B, therefore your final operation after implementing the three basic methods should be:

AT*BT*(B*BT)'

Note: Don't forget the basic rules of operator precedence :)

2009-06-16 14:07
by Suvesh Pratapa


1

[edited] As Suvesh pointed out, i was completely wrong before. however, numpy can still easily do the procedure he gives in his post:

A = numpy.matrix(numpy.random.random((18, 1))) # as noted by others, your dimensions are off
B = numpy.matrix(numpy.random.random((25, 18)))
C = A.T * B.T * (B * B.T).I
2009-06-16 14:34
by Autoplectic
The / operator in Python is defined as standard matrix division for square matrices, i.e., A*inv(B). In her example, she's trying to achieve right division for any size matrices. So your code would not work - Suvesh Pratapa 2009-06-16 14:38
+1 for the change, and translating my math into Python code. :)

I don't know Python, but in any case, I wanted her to work that out, so I didn't bother to post code along with my math - Suvesh Pratapa 2009-06-16 15:20

This at least gives me a result matrix c with the right dimensions, but the values in that matrix do not match the values in Matlab- any ideas - EmilyS 2009-06-16 15:26
did you make sure that the A and B you use in numpy are the same A and B you used in matlab - Autoplectic 2009-06-16 16:02
You should NOT do this. The so called normal equations should not be implemented directly, but should use the least square approximation, which is lstsq function in numpy.linalg - David Cournapeau 2009-06-16 16:04
@Autoplectic: yes. @David: can you please explain how to use this function to replace this line in Python? I tried to use it after I read the help for lstsq, but I am getting getting values in the matrix that are different from the ones in Matlab, though they are the same dimensions - EmilyS 2009-06-16 17:28
There is rarely any good reason to invert a matrix. See, for instance, this excellent blog post. Also, there's rarely any good reason to use numpy.matrix instead of numpy.ndarray. The docs give several compelling reasons to favor ndarrays over matrixes - Praveen 2018-03-07 18:26


1

You can also approach this using the pseudo-inverse of B then post multiplying that result with A. Try using numpy.linalg.pinv then combine this with matrix multiplication via numpy.dot:

c = numpy.dot(a, numpy.linalg.pinv(b))
2018-03-07 16:20
by rayryeng
Ads