CODEWITHSUNDEEP
pythonarraysnumpymatrix
Josh Pinto
Josh Pinto

Reputation: 1573

Overcoming dimension errors with matrix multiplication?

I'd like to multiply two matrices (to be clear I want it in the linear algebra meaning of the word, i.e., if you want further clarification of what I mean see this Wikipedia article) I've tried A * B (which I got by following the MATLAB to Numpy guide at Scipy Wiki) but when A had dimensions of (N-1)x(N-1) and B had dimensions of (N-1)x(N+1) it gave out the ValueError (where N=100):

ValueError: operands could not be broadcast together with shapes (99,99) (99,101)

I also tried np.dot(A,B) in case I was unintentionally using an array (ps. this could be the case, as I do have a few np.asarray commands for the vectors I built A and B from. How do I convert them to matrices?) and not a matrix and got the same result.

This is my code in full (granted A and B are replaced by np.diag(xasub) and np.sin..., respectively and this is on the d2Tsub line)

import numpy as np
import scipy.special as sp

N = 100
n = range(0,N+1)
xmin = -10
xmax = 5

pi = np.pi
x = np.cos(np.multiply(pi / float(N), n))

xa = np.asarray(x)
xasub = xa[1:N]

xc = (1+xa)*(xmax-xmin) / float(2)
xcsub = xc[1:N]

na = np.asarray(n)
nd = np.transpose(na)

T = np.cos(np.outer(np.arccos(xa),nd))
Tsub  = T[1:-1]

from numpy import linalg as LA

dT = np.sin(np.outer(np.arccos(xa),nd))*na
dTsub = dT[1:-1]

d2Tsub = np.diag(np.power((xasub*xasub-1),-1)) * (Tsub*nd - np.diag(xasub) * np.sin(np.outer(np.arccos(xasub),nd))) * nd

Upvotes: 0

Views: 301

Answers (2)

hpaulj
hpaulj

Reputation: 231385

I've replicated your MATLAB code with:

N = 100.
n = np.arange(N+1).reshape(1,-1)
x = np.cos(np.pi*n/N).T
xsub = x[1:-1,:]
T = np.cos(np.arccos(x)*n)  # broadcasted outer product
Tsub = T[1:-1,:]
dT = np.dot(np.sin(np.arccos(x)*n), np.diag(n.flat))
# np.diag needs 1d input to generate 2d array
# this is a case of matrix multiplication, the np.dot for array
# * for np.matrix

continuing

temp = np.dot(np.diag(xsub.flat), np.sin(np.arccos(xsub)*n))
temp = np.dot(Tsub, np.diag(n.flat)) - temp
temp = np.dot(temp, np.diag(n.flat))
d2Tsub = np.dot(np.diag((1/np.sqrt(1-xsub**2)).flat), temp))

messy, but runs and generates matching value. Use of np.matrix versions of these arrays might simplify the appearance, as would cleaning up the repeated np.diag(...flat).

for the rest of d2T, this vstack is a pretty literal translation.

d2T1 = ((-1)**n) * (n**2-1) * (n**2)/3
d2T2 = (n**2-1) * (n**2)/3
d2T = np.vstack([d2T1, d2Tsub, d2T2])

last step is a matrix divide. MATLAB has 2 versions:

D2 = np.linalg.solve(T, d2T)        # matlab T\d2T
D2 = np.linalg.solve(T.T, d2T.T).T  # matlab d2T/T

I worked out this last equivalence from the Octave docs for Y/X and X\Y .

This is a literal translation. I'm not trying to understand what is going on, or to make it idiomatic numpy.

A partial simplification. No need for matrix products to use diagonal matrices. Simple (broadcasted) element multiplication is enough.

n = np.arange(N+1)
x = np.cos(np.pi*n/N)[:,None]
Xn = np.arccos(x)*n
T = np.cos(Xn)
dT = np.sin(Xn)*n

xsub = x[1:-1,:]
d2Tsub = ((1-xsub**2)**(-.5)) * (T[1:-1,:] * n - xsub * np.sin(Xn[1:-1,:]))
d2T2 = (n**2-1) * (n**2)/3
d2T1 = ((-1)**n) * d2T2

You could probably replace a lot of the diag(n)* in matlab with n.* For example

dT = n.*sin(acos(x)*n);  # equals
dT = sin(acos(x)*n)*diag(n);

However, I tested this in Octave which now has broadcasting. I'm not sure where MATLAB stands in that regard.

Upvotes: 1

jedwards
jedwards

Reputation: 30210

In fact, you have all arrays and no matrices:

xasub : <type 'numpy.ndarray'>
xc    : <type 'numpy.ndarray'>
xa    : <type 'numpy.ndarray'>
xmin  : <type 'int'>
LA    : <type 'module'>
na    : <type 'numpy.ndarray'>
nd    : <type 'numpy.ndarray'>
Tsub  : <type 'numpy.ndarray'>
np    : <type 'module'>
pi    : <type 'float'>
dTsub : <type 'numpy.ndarray'>
N     : <type 'int'>
T     : <type 'numpy.ndarray'>
sp    : <type 'module'>
xcsub : <type 'numpy.ndarray'>
n     : <type 'list'>
xmax  : <type 'int'>
x     : <type 'numpy.ndarray'>
dT    : <type 'numpy.ndarray'>

So you'll need to modify the last line to use np.dot() accordingly, for example (with horrible variable names):

#d2Tsub = np.diag(np.power((xasub*xasub-1),-1)) * (Tsub*nd - np.diag(xasub) * np.sin(np.outer(np.arccos(xasub),nd)))*nd * nd
d2T_1  = np.diag(np.power((xasub*xasub-1),-1))
d2T_2  = np.dot( np.diag(xasub), np.sin(np.outer(np.arccos(xasub),nd)) )

d2T_3  = np.dot(d2T_1, (Tsub*nd - d2T_2))

d2Tsub = np.dot(d2T_3, nd)

print d2Tsub.shape  # (99L,)

Of course, you could use the a.dot(b) syntax instead of np.dot(a,b) and save a tiny bit space.

Upvotes: 0

Related Questions