# inverting a matrix using the LU decomposition approach

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.

• 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...) – Ian Jan 24 '18 at 15:21
• 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. – Yves Daoust Jan 24 '18 at 15:22
• @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. – Dalek Jan 24 '18 at 15:26
• 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). – Ian Jan 24 '18 at 15:28
• 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.) – Ian Jan 24 '18 at 15:33