 <p>You could use a sparse block matrix A which stores the (5, 2) entries of T_Arm on its diagonal, and solve AX = b where b is the vector composed of stacked entries of <code>Erg</code>. Then solve the system with scipy.sparse.linalg.lsqr(A, b).</p> <p>To construct A and b I use n=3 for visualisation purposes:</p> <pre><code>import numpy as np import scipy from scipy.sparse import bsr_matrix n = 3 col = np.hstack(5 * [np.arange(10 * n / 5).reshape(n, 2)]).flatten() array([ 0., 1., 0., 1., 0., 1., 0., 1., 0., 1., 2., 3., 2., 3., 2., 3., 2., 3., 2., 3., 4., 5., 4., 5., 4., 5., 4., 5., 4., 5.]) row = np.tile(np.arange(10 * n / 2), (2, 1)).T.flatten() array([ 0., 0., 1., 1., 2., 2., 3., 3., 4., 4., 5., 5., 6., 6., 7., 7., 8., 8., 9., 9., 10., 10., 11., 11., 12., 12., 13., 13., 14., 14.]) A = bsr_matrix((T_Arm[:n].flatten(), (row, col)), shape=(5 * n, 2 * n)) A.toarray() array([[ 0, 1, 0, 0, 0, 0], [ 2, 3, 0, 0, 0, 0], [ 4, 5, 0, 0, 0, 0], [ 6, 7, 0, 0, 0, 0], [ 8, 9, 0, 0, 0, 0], [ 0, 0, 10, 11, 0, 0], [ 0, 0, 12, 13, 0, 0], [ 0, 0, 14, 15, 0, 0], [ 0, 0, 16, 17, 0, 0], [ 0, 0, 18, 19, 0, 0], [ 0, 0, 0, 0, 20, 21], [ 0, 0, 0, 0, 22, 23], [ 0, 0, 0, 0, 24, 25], [ 0, 0, 0, 0, 26, 27], [ 0, 0, 0, 0, 28, 29]], dtype=int64) b = Erg[:n].flatten() </code></pre> <p>And then </p> <pre><code>scipy.sparse.linalg.lsqr(A, b)[0] array([ 5.00000000e-01, -1.39548109e-14, 5.00000000e-01, 8.71088538e-16, 5.00000000e-01, 2.35398726e-15]) </code></pre> <p>EDIT: A is not as huge in memory as it seems: more on block sparse matrices <a href="https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.bsr_matrix.html" rel="nofollow noreferrer">here</a>.</p>