Ada Xu
Ada Xu

Reputation: 983

Multiplying Block Matrices in Numpy

Hi Everyone I am python newbie I have to implement lasso L1 regression for a class assignment. This involves solving a quadratic equation involving block matrices.

minimize x^t * H * x  + f^t * x 
where x > 0

Where H is a 2 X 2 block matrix with each element being a k dimensional matrix and x and f being a 2 X 1 vectors each element being a k dimension vector.

I was thinking of using ndarrays.

Such that :

  np.shape(H) = (2, 2, k, k)
  np.shape(x) = (2, k)

But I figured out that np.dot(X, H) doesn't work here. Is there an easy way to solve this problem? Thanks in advance.

Upvotes: 5

Views: 4006

Answers (2)

hpaulj
hpaulj

Reputation: 231530

np.einsum gives good control over which axes are summed.

np.einsum('ijkl,jk',H,x)

is one possible (generalized) dot product, (2,4) (first and last dim of H)

np.einsum('ijkl,jl',H,x)

is another. You need to be explicit - which dimensions of x go with which of H.

Upvotes: 1

alko
alko

Reputation: 48347

First of all, I am convinced that converting to matrices will lead to more efficient computations. Stating that, if you consider your 2k x 2k matrix being a 2 x 2 matrix, then you operate in a tensor product of vector spaces, and have to use tensordot instead of dot.

Let give it a try, with k=5 for example:

>>> import numpy as np
>>> k = 5

Define our matrix a and vector x

>>> a = np.arange(1.*2*2*k*k).reshape(2,2,k,k)
>>> x = np.arange(1.*2*k).reshape(2,k)
>>> x
array([[ 0.,  1.,  2.,  3.,  4.],
       [ 5.,  6.,  7.,  8.,  9.]])

now we can multipy our tensors. Be sure to choose right axes, I didn't tested following formula explicetely, and there might be an error

>>> result = np.tensordot(a,x,([1,3],[0,1]))
>>> result
array([[  985.,  1210.,  1435.,  1660.,  1885.],
       [ 3235.,  3460.,  3685.,  3910.,  4135.]])
>>> np.shape(result)
(2, 5)

Upvotes: 1

Related Questions