2
$\begingroup$

I have written the Matrix class in cython for the matrix inversion and some other linear algebra operations. I tried to use the LU decomposition and forward, backward substitutions, in order to measure the inverse of my matrix. The speed of code is good, but somehow the elements of the inverted matrix are mixed up and I could not figure out what I have done wrong. I will appreciate for any help to find the mistake I made in my code.

matrix.pyx

cdef class Matrix:    
     def __cinit__(self, size_t rows=0, size_t columns=0, bint Identity=False, bint ones=False):
         self._rows=rows
         self._columns=columns
         self.matrix=new vector[double]()
         self.matrix.resize(rows*columns)

         if Identity:
            self._IdentityMatrix()

         if ones:
            self._fillWithOnes()

     def __dealloc__(self):
         del self.matrix

     property rows:
        def __get__(self):
            return self._rows
        def __set__(self, size_t x):
            self._rows = x    
     property columns:
        def __get__(self):
            return self._columns
        def __set__(self, size_t y):
            self._columns = y    

     cpdef double getVal(self, size_t r, size_t c):
           return self.matrix[0][r*self._columns+c]
     cpdef void setVal(self, size_t r, size_t c, double v): 
           self.matrix[0][r*self._columns+c] = v
     @cython.boundscheck(False)
     @cython.wraparound(False)
     cdef void _fillWithOnes(self):
          fill(self.matrix.begin(),self.matrix.end(),1.)

     cdef void _IdentityMatrix(self):
          cdef size_t i 
          if (self._rows!=self._columns):
             raise Exception('In order to generate identity matrix, the number of rows and columns must be equal')
          else:
             for i from 0 <= i <self._columns:
                 self.setVal(i,i,1.0)

     @cython.boundscheck(False)
     @cython.wraparound(False)
     cpdef Matrix Inv(self):
           #https://github.com/Kamp9/Bachelorprojekt/blob/93067b9bcf9724ca129afb7735b05111f7871278/solve.py
           cdef Matrix B      = Matrix(self._rows,self._columns, Identity=True)
           cdef MatrixList LU = lu_decomposition(self)
           cdef Matrix L      = LU.get(0)
           cdef Matrix U      = LU.get(1)
           cdef int j,i
           cdef Matrix Z, X
           cdef Matrix R=Matrix(self._rows,1)
           cdef Matrix A_inverse = Matrix(self._rows,self._columns)
           for j from 0 <= j < self._columns:
               for i from 0 <= i < self._rows:
                   R.setVal(i,0,B.getVal(i,j))
               Z = ForwardSubstitution(R, L)
               X = BackSubstitution(U, Z)
               for i from 0 <= i < self._rows:
                   A_inverse.setVal(i,j,X.getVal(i,0))
           return A_inverse

@cython.boundscheck(False)
@cython.wraparound(False)    
cdef Matrix pivotize(Matrix m):
      """Creates the pivoting matrix for m."""
      cdef int n = m.rows
      cdef int i, j, row, index
      # Create an identity matrix, with floating point values 
      cdef Matrix ID = Matrix(n,n, Identity= True)
      cdef Matrix tmp = Matrix(1,n)
      cdef double maxVal, mmin, mmax
      cdef vector[double].iterator vin=min_element(m.matrix.begin(), m.matrix.end())
      index = vin - m.matrix.begin()
      mmin  = m.matrix.at(index)
      maxVal = mmin
      for j from 0 <= j < n:
          mmax=maxVal
          for i from j <= i < n:
              if (fabs(m.getVal(i,j))>mmax):
                 row = i
                 mmax = fabs(m.getVal(i,j))
          if j != row:
             copy(ID.matrix.begin() + j*ID._columns , ID.matrix.begin() + (j+1)*ID._columns, tmp.matrix.begin())
             copy(ID.matrix.begin() + row*ID._columns, ID.matrix.begin() +(row+1)*ID._columns, ID.matrix.begin()+j*ID._columns)
             copy(tmp.matrix.begin(), tmp.matrix.end(),ID.matrix.begin()+row*ID._columns)
      return ID

cdef class MatrixList:
     def __cinit__(self):
         self.inner = []

     cdef void append(self, Matrix a):
          self.inner.append(a)

     cdef Matrix get(self, int i):
          return <Matrix> self.inner[i]

     def __len__(self):
         return len(self.inner)

@cython.boundscheck(False)
@cython.wraparound(False)    
cdef Matrix ForwardSubstitution(Matrix B, Matrix L):
     # Forward solve Ly = b
     cdef int n    = B._rows
     cdef Matrix Y = Matrix(n,1)
     cdef int i,j
     for i from 0 <= i < n:
         Y.setVal(i,0,B.getVal(i,0))
         for j from 0 <= j < i:
             Y.setVal(i,0,Y.getVal(i,0)-L.getVal(i,j)*Y.getVal(j,0))
         Y.setVal(i,0,Y.getVal(i,0)/L.getVal(i,i))
     return Y

@cython.boundscheck(False)
@cython.wraparound(False)    
cdef Matrix BackSubstitution(Matrix U, Matrix z): 
      #Backward solve Ux = y
      #https://github.com/Kyezil/ALN/blob/7ff209651d2ef577af845c101ca96e39d3cb38f4/PracticaC%2B%2B/src/LU_decomposition.cc
      cdef int m=U._rows
      cdef int n=U._columns
      cdef Matrix X =Matrix(n,1)
      cdef int i,j
      cdef double s
      X.setVal(n-1,0,z.getVal(n-1,0)/U.getVal(n-1,n-1))
      for i from n-1 > i >= 0:
          s=0.0
          for j from i+1 <= j < n:
              s+= X.getVal(j,0)*U.getVal(i,j)
          X.setVal(i,0, (z.getVal(i,0)- s)/U.getVal(i,i))  
      return X

@cython.boundscheck(False)
@cython.wraparound(False)    
cdef MatrixList lu_decomposition(Matrix A):
    #https://www.quantstart.com/articles/LU-Decomposition-in-Python-and-NumPy
    """Performs an LU Decomposition of A (which must be square)                                                                                                                                                                                        
    into PA = LU. The function returns P, L and U."""
    cdef int n = A._rows
    cdef int i,j,k
    # Create zero matrices for L and U                                                                                                                                                                                                                 
    cdef Matrix L      = Matrix(n,n)
    cdef Matrix U      = Matrix(n,n)
    cdef double s1, s2
    # Create the pivot matrix P and the multipled matrix PA                                                                                                                                                                                            
    cdef Matrix P      = pivotize(A)
    cdef Matrix PA     = MatrixMultiplication(P, A)
    cdef MatrixList LU = MatrixList()
    for j from 0 <= j < n:
        L.setVal(j,j,1.)
        for i from 0 <= i < (j+1):
            s1=0.0
            for k from 0 <= k < i:
                s1 += U.getVal(k,j)*L.getVal(i,k)
            U.setVal(i,j, PA.getVal(i,j)-s1)
        for i from j <= i < n:        
            s2=0.0
            for k from 0 <= k < j:
                s2 += U.getVal(k,j)*L.getVal(i,k)
            L.setVal(i,j, (PA.getVal(i,j)-s2)/U.getVal(j,j))
    LU.append(<Matrix>L)
    LU.append(<Matrix>U)
    LU.append(<Matrix>P)
    return LU        


@cython.boundscheck(False)
@cython.wraparound(False)    
cpdef Matrix MatrixMultiplication(Matrix matrix1, Matrix matrix2):
      cdef size_t r1 = matrix1._rows
      cdef size_t c1 = matrix1._columns
      cdef size_t r2 = matrix2._rows
      cdef size_t c2 = matrix2._columns
      cdef Matrix multiplication = Matrix(r1, c2)
      cdef size_t row, colIndex, col
      cdef double temp
      if (c1 == r2):
         for row from 0 <= row < r1:
             for colIndex from 0 <= colIndex <c2:
                 # Iterate rows index of second matrix
                 for col from 0 <= col <r2:
                     temp = multiplication.getVal(row,colIndex)   
                     multiplication.setVal(row, colIndex, temp + matrix1.getVal(row,col) * matrix2.getVal(col,colIndex))
      else:
         raise Exception("Multiplication not Possible; required columns of first matrix equals to the rows of second matrix.")
      return multiplication

matrix.pxd

from libcpp.vector cimport vector
cdef class MatrixList:
     cdef list inner
     cdef void append(self, Matrix a)
     cdef Matrix get(self, int i)

cdef class Matrix:
     cdef vector[double] *matrix   
     cdef size_t _rows
     cdef size_t _columns
     cdef bint Identity
     cdef bint ones

     cpdef double getVal(self, size_t r, size_t c)
     cpdef void setVal(self, size_t r, size_t c, double v)
     cpdef Matrix transpose(self)
     cdef void _IdentityMatrix(self)
     cdef void _fillWithOnes(self)
     cpdef Matrix Inv(self)

cdef Matrix pivotize(Matrix m)
cdef Matrix ForwardSubstitution(Matrix B, Matrix L)
cdef Matrix BackSubstitution(Matrix U, Matrix z)
cdef MatrixList lu_decomposition(Matrix A)
cpdef Matrix MatrixMultiplication(Matrix matrix1, Matrix matrix2)

I checked my LU decomposition function and works correctly, I guess my main mistake must occurred in Inv() function of Matrix class where I fill the inverted matrix. But I am also not 100% sure about the BackSubstitution and ForwardSubstitution functions. These are my three suspicious points in the code.

>>>import numpy
>>>from matrix import Matrix
>>>from numpy.linalg import inv
>>>import timeit
>>>import numpy as np
>>>r=numpy.random.random((20, 20))

>>>d=Matrix(r.shape[0],r.shape[1])
>>> for i in range(d.rows):
        for j in range(d.columns):
            d.setVal(i,j,r[i,j])
            print d.getVal(i,j),r[i,j]

>>> start = timeit.default_timer()
>>> x=d.Inv()
>>> stop = timeit.default_timer()
>>> print "LU decomposition:", stop - start 
LU decomposition: 0.000540018081665
>>> start = timeit.default_timer()
>>> a=inv(r)
>>> stop = timeit.default_timer()
>>> print "numpy inverse :", stop - start 
numpy inverse : 0.00105309486389

I did a test for different size matrices and for matrix with row=column=3 the first column swapped with the second one.

$\endgroup$
  • $\begingroup$ Directly obtaining the actual inverse from the LU decomposition is probably not as simple as you think. Indeed it appears to me that the main way to do this is to solve $LUx=b$ for $b=e_i$ for $i=1,2,\dots,n$, which is quite slow (a bit harder than the LU decomposition was to compute in the first place). But you rarely need the actual inverse, in practice if you have the LU decomposition then you can solve $Ax=b$ by solving $Ly=b$ and $Ux=y$, each by back-substitution. (I will also add that you are reinventing the wheel here, but presumably there is some motivation for that...) $\endgroup$ – Ian Jan 24 '18 at 15:21
  • $\begingroup$ SOrry, this is not an answer, more a criticism: it is a pity to use so many matrices when the LU method has been designed to be performed in-place, with only constant extra space. $\endgroup$ – Yves Daoust Jan 24 '18 at 15:22
  • $\begingroup$ @Ian I used exactly the approach you have suggested and for $b$ I used the identity matrix columns. I think in this process, I mixed up things. Well, I needed some better speed over numpy which I reached plus, in my class I do other matrix operations. $\endgroup$ – Dalek Jan 24 '18 at 15:26
  • $\begingroup$ I'm not sure how, then: you should just have $(A^{-1})_i=U^{-1} L^{-1} e_i$. But your code is a little bit hard to follow (because there is a bit of mixing of high-level logic and low-level logic, if you understand my meaning). $\endgroup$ – Ian Jan 24 '18 at 15:28
  • $\begingroup$ I would also be surprised if this was actually faster than numpy, if numpy is installed with optimization in mind (linked to locally compiled BLAS etc.) $\endgroup$ – Ian Jan 24 '18 at 15:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.