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.