3

I was wondering if there is an easy way to assemble matrices along the diagonal in python, adding values if they overlap. Here is a handy diagram I stole from the Matlab forum: https://i.stack.imgur.com/ZunrZ.jpg

Currently, I am trying to implement it to work with a set number of 2x2 matrices, but the end goal is to make the code assemble an arbitrary amount of arbitrary sized matrices(all the same size ofc., max 4x4).

Iguananaut
  • 21,810
  • 5
  • 50
  • 63
Preslav Aleksandrov
  • 31
  • 1
  • 3
  • How do you want the overlap to be if they are 3x3 submatrices? 1 element overlap, or 4 elements overlap? – wim May 07 '18 at 22:05
  • I'm not quite sure... in the model, I'm creating I planned to only have 1 overlap. But it might turn out to do the thing I'm doing I might need more complex scenarios. Nevertheless, for now, just one corner overlap is what I'm looking for. – Preslav Aleksandrov May 07 '18 at 22:10
  • This looks a lot like the assembly of a global stiffness matrix from element ones for finite element modeling. `scipy.sparse` does this when converting from `coo` to `csr` format matrices, in much the same as a MATLAB does with its sparse matrices. – hpaulj May 07 '18 at 23:00
  • See references to finite element construction on this page: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html#scipy.sparse.coo_matrix – hpaulj May 07 '18 at 23:05
  • 1
    This looks related: https://stackoverflow.com/questions/21139822/build-scipy-matrix-with-diagonal-overlays – Warren Weckesser May 07 '18 at 23:30
  • 1
    And [Compiling n submatrices into an NxN matrix in numpy](https://stackoverflow.com/questions/37039923/compiling-n-submatrices-into-an-nxn-matrix-in-numpy) – hpaulj May 08 '18 at 01:35

2 Answers2

3

I'm not sure how to vectorize this, but you can do it fairly directly with setitem:

k = k1 = np.array([[1,2],[3,4]])   # etc
ks = [k1, k2, k3, k4]
[n] = set(k.shape)
N = len(ks)
A = np.zeros((N+1, N+1))
for i, k in enumerate(ks):
    A[i:i+n, i:i+n] += k
wim
  • 338,267
  • 99
  • 616
  • 750
  • Thank you so much!!! Quick follow up does np have a way of storing/assigning the location of each submatrix within itself so that I can assemble them out of order or in a different direction? i.e instead of using the index of the matrix for locating to use a value stored within the matrix. – Preslav Aleksandrov May 07 '18 at 22:20
  • No, it does not - you would have to waste a whole row and column for this to store it in the matrix itself. Better would be to use a list of pairs (ndarray, idx). – wim May 07 '18 at 22:29
  • Thanks again, that’s perfect for what I’m doing. – Preslav Aleksandrov May 07 '18 at 22:32
1

First an example with two 2x2 matrices on a 5x5 big matrix:

import numpy as np
M = np.zeros((5,5))
M1 = np.matrix([[1,2],[3,4]])
M2 = np.matrix([[1,2],[3,4]])
M[:2,:2] += M1
M[1:3, 1:3] += M2
M

yields

array([[ 1.,  2.,  0.,  0.,  0.],
   [ 3.,  5.,  2.,  0.,  0.],
   [ 0.,  3.,  4.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]])

In general, if you have a KxK matrix of zeros M, and K-r+1 rxr matrices in some indexable, you can do

for i in range(K-r+1): 
    M[i+r:i+r] += matrices[i]
fabikw
  • 109
  • 4