Solve a multitude of linear least square system efficiently

Question

I have to find the best solution for >10^7 equation systems with 5 equations in 2 variables each (5 measurements to find 2 parameters with the least amount of error in a long series). The following code (normally used to do curve fitting) does what I want:

#Create_example_Data
n = 100
T_Arm = np.arange(10*n).reshape(-1, 5, 2)
Erg = np.arange(5*n).reshape(-1, 5)
m = np.zeros(n)
c = np.zeros(n)
#Run
for counter in xrange(n):
     m[counter], c[counter] = np.linalg.lstsq(T_Arm[counter, :, :], 
                                              Erg[counter, :])[0]

Unfortunately it is too slow. Is there any way how to speed this code up significantly? I tried to vectorise it, but I did not succeed. Using the last solution as a initial guess might be a good idea as well. Using scipy.optimize.leastsq did not speed it up as well.


Show source
| numpy   | python   | performance   | linear-algebra   2016-11-09 13:11 1 Answers

Answers ( 1 )

  1. 2016-11-09 16:11

    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 Erg. Then solve the system with scipy.sparse.linalg.lsqr(A, b).

    To construct A and b I use n=3 for visualisation purposes:

    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()
    

    And then

    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])
    

    EDIT: A is not as huge in memory as it seems: more on block sparse matrices here.

◀ Go back