[ACCEPTED]-Array division- translating from MATLAB to Python-linear-algebra

Accepted answer
Score: 10

The line

c = a.' / b

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

c = numpy.linalg.lstsq(b.T, a.T)[0].T
Score: 9

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

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

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

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

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

Score: 5

In Matlab, A.' means transposing the A matrix. So 11 mathematically, what is achieved in the 10 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 9 (You want to find C as A/B)

RIGHT DIVISION 8 (/) is as follows:

C*(B*BT)=A*BT

You then isolate 7 C by inverting (B*BT)

i.e.,

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

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

  • Matrix multiplication
  • Matrix transpose
  • Matrix inverse

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

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

AT*BT*(B*BT)'

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

Score: 2

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

c = numpy.dot(a, numpy.linalg.pinv(b))
Score: 1

[edited] As Suvesh pointed out, i was completely 3 wrong before. however, numpy can still 2 easily do the procedure he gives in his 1 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

More Related questions