使用 LU 分解求逆矩阵
Inverting a matrix using the LU decomposition
我在 cython 中编写了以下 Matrix class 用于矩阵求逆和其他一些线性代数运算。我尝试使用 LU 分解 来计算矩阵的逆。代码速度很好。我试图在 cython 中实现 this code。我已经检查了我的代码的每一行并与给定的代码进行了几次比较，但我仍然 return 错误的答案。
matrix.pyx
from libcpp.vector cimport vector
import cython
cimport cython
import numpy as np
cimport numpy as np
import ctypes
from libc.math cimport log, exp, pow, fabs
from libc.stdint cimport *
from libcpp.string cimport string
from libc.stdio cimport *
from libcpp cimport bool
cdef extern from "<iterator>" namespace "std" nogil:
cdef cppclass iterator[Category, T, Distance, Pointer, Reference]:
pass
cdef cppclass output_iterator_tag:
pass
cdef cppclass input_iterator_tag:
pass
cdef cppclass forward_iterator_tag(input_iterator_tag):
pass
cdef extern from "<algorithm>" namespace "std" nogil:
void fill [ForwardIterator, T](ForwardIterator, ForwardIterator, T& )
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):
cdef Matrix A_inverse = Matrix(self._rows,self._columns)
cdef MatrixList LU = ludcmp(self)
cdef Matrix A = LU.get(0)
cdef Matrix indx = LU.get(1)
cdef Matrix d = LU.get(2)
cdef double det = d.getVal(0,0)
cdef int i, j
cdef np.ndarray[np.float64_t, ndim=2] L = np.zeros((self._rows,self._columns),dtype=np.float64)
cdef np.ndarray[np.float64_t, ndim=2] U = np.zeros((self._rows,self._columns),dtype=np.float64)
cdef Matrix col = Matrix(self._rows,1)
for i from 0 <= i < self._rows:
for j from 0 <= j < self._columns:
if (j>i):
U[i,j]=A.getVal(i,j)
L[i,j]=0
elif (j<i):
U[i,j]=0
L[i,j]=A.getVal(i,j)
else:
U[i,j]=A.getVal(i,j)
L[i,j]=1
print "product of a lower triangular matrix L and an upper triangular matrix U:", np.dot(L, U)
for i from 0 <= i < self._rows:
det*= A.getVal(i,i)
for j from 0 <= j < self._columns:
if (i==j):
col.setVal(j,0,1)
col=lubksb(A, indx, col)
for j from 0 <= j < self._columns:
A_inverse.setVal(j,i,col.getVal(j,0))
print "determinant of matrix %.4f"%(det)
return A_inverse
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 lubksb(Matrix a, Matrix indx, Matrix b):
cdef int n = a.rows
cdef int i, ip, j
cdef int ii = 0
cdef double su
for i from 0 <= i < n:
ip = <int>indx.getVal(i,0)
su = b.getVal(ip,0)
b.setVal(ip,0, b.getVal(i,0))
if (ii):
for j from ii <= j < (i-1):
su -= a.getVal(i,j) * b.getVal(j,0)
elif (su):
ii = i
b.setVal(i, 0, su)
for i from n > i >= 0:
su = b.getVal(i,0)
for j from (i+1) <= j < n:
su -= a.getVal(i,j) * b.getVal(j,0)
b.setVal(i, 0, su/a.getVal(i,i))
return b
@cython.boundscheck(False)
@cython.wraparound(False)
cdef MatrixList ludcmp(Matrix a):
#Given a matrix a_{nxn}, this routine replaces it by the LU decomposition of a row-wise permutation of itself.
cdef MatrixList LU = MatrixList()
cdef int n = a.rows
cdef int i, j, k, imax
cdef double big, dum, su, temp
cdef Matrix vv = Matrix(n,1)
cdef Matrix indx = Matrix(n,1) #an output vector that records the row permutation effected by the partial pivoting
cdef Matrix d = Matrix(1,1, ones= True) #an output as +1 or -1 depending on whether the number of row interchanges was even or odd, respectively
cdef double TINY = 1.1e-16
for i from 0 <= i < n:
big = 0.0
for j from 0 <= j < n:
temp=fabs(a.getVal(i,j))
if (temp > big):
big=temp
if (big ==0.0):
raise Exception("ERROR! ludcmp: Singular matrix\n")
vv.setVal(i,0,1.0/big)
for j from 0 <= j < n:
for i from 0 <= i < j:
su = a.getVal(i,j)
for k from 0 <= k < i:
su -= a.getVal(i,k)*a.getVal(k,j)
a.setVal(i,j,su)
big=0.0
for i from j<= i< n:
su = a.getVal(i,j)
for k from 0 <= k < j:
su -= a.getVal(i,k)*a.getVal(k,j)
a.setVal(i, j, su)
dum=vv.getVal(i,0)*fabs(su )
if (dum >= big):
big=dum
imax=i
if (j != imax):
for k from 0 <= k < n:
dum = a.getVal(imax,k)
a.setVal(imax, k, a.getVal(j,k))
a.setVal(j,k, dum)
d.setVal(0, 0, -d.getVal(0,0))
vv.setVal(imax, 0, vv.getVal(j, 0))
indx.setVal(j, 0, imax)
if (a.getVal(j,j) == 0.0):
a.setVal(j,j, TINY)
if (j != (n-1)):
dum=1.0/a.getVal(j,j)
for i from (j+1)<= i <n:
a.setVal(i,j, a.getVal(i,j)*dum)
LU.append(<Matrix>a)
LU.append(<Matrix>indx)
LU.append(<Matrix>d)
return LU
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 lubksb(Matrix a, Matrix indx, Matrix b)
cdef MatrixList ludcmp(Matrix a)
任何帮助查找错误的帮助将不胜感激。
示例：
import numpy
from matrix import Matrix
from numpy.linalg import inv
import timeit
import numpy as np
r=numpy.random.random((100, 100))
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])
start = timeit.default_timer()
x=d.Inv()
stop = timeit.default_timer()
print "LU decomposition:", stop - start
我发现我在 lubksb 函数中犯了一些非常小的错误，通过修复这些错误我得到了正确的答案。这是固定代码：
@cython.boundscheck(False)
@cython.wraparound(False)
cdef Matrix lubksb(Matrix a, Matrix indx, Matrix b):
cdef int n = a.rows
cdef int i, ip, j
cdef int ii = 0
cdef double su
for i from 0 <= i < n:
ip = <int>indx.getVal(i,0)
su = b.getVal(ip,0)
b.setVal(ip,0, b.getVal(i,0))
if (ii>=0):
for j from ii <= j <= (i-1):
su -= a.getVal(i,j) * b.getVal(j,0)
elif (su):
ii = i
b.setVal(i, 0, su)
for i from n >= i >= 0:
su = b.getVal(i,0)
for j from (i+1) <= j < n:
su -= a.getVal(i,j) * b.getVal(j,0)
if (a.getVal(i,i)==0.0):
a.setVal(i,i, 1.1e-16)
b.setVal(i, 0, su/a.getVal(i,i))
return b
我在 cython 中编写了以下 Matrix class 用于矩阵求逆和其他一些线性代数运算。我尝试使用 LU 分解 来计算矩阵的逆。代码速度很好。我试图在 cython 中实现 this code。我已经检查了我的代码的每一行并与给定的代码进行了几次比较，但我仍然 return 错误的答案。
matrix.pyx
from libcpp.vector cimport vector
import cython
cimport cython
import numpy as np
cimport numpy as np
import ctypes
from libc.math cimport log, exp, pow, fabs
from libc.stdint cimport *
from libcpp.string cimport string
from libc.stdio cimport *
from libcpp cimport bool
cdef extern from "<iterator>" namespace "std" nogil:
cdef cppclass iterator[Category, T, Distance, Pointer, Reference]:
pass
cdef cppclass output_iterator_tag:
pass
cdef cppclass input_iterator_tag:
pass
cdef cppclass forward_iterator_tag(input_iterator_tag):
pass
cdef extern from "<algorithm>" namespace "std" nogil:
void fill [ForwardIterator, T](ForwardIterator, ForwardIterator, T& )
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):
cdef Matrix A_inverse = Matrix(self._rows,self._columns)
cdef MatrixList LU = ludcmp(self)
cdef Matrix A = LU.get(0)
cdef Matrix indx = LU.get(1)
cdef Matrix d = LU.get(2)
cdef double det = d.getVal(0,0)
cdef int i, j
cdef np.ndarray[np.float64_t, ndim=2] L = np.zeros((self._rows,self._columns),dtype=np.float64)
cdef np.ndarray[np.float64_t, ndim=2] U = np.zeros((self._rows,self._columns),dtype=np.float64)
cdef Matrix col = Matrix(self._rows,1)
for i from 0 <= i < self._rows:
for j from 0 <= j < self._columns:
if (j>i):
U[i,j]=A.getVal(i,j)
L[i,j]=0
elif (j<i):
U[i,j]=0
L[i,j]=A.getVal(i,j)
else:
U[i,j]=A.getVal(i,j)
L[i,j]=1
print "product of a lower triangular matrix L and an upper triangular matrix U:", np.dot(L, U)
for i from 0 <= i < self._rows:
det*= A.getVal(i,i)
for j from 0 <= j < self._columns:
if (i==j):
col.setVal(j,0,1)
col=lubksb(A, indx, col)
for j from 0 <= j < self._columns:
A_inverse.setVal(j,i,col.getVal(j,0))
print "determinant of matrix %.4f"%(det)
return A_inverse
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 lubksb(Matrix a, Matrix indx, Matrix b):
cdef int n = a.rows
cdef int i, ip, j
cdef int ii = 0
cdef double su
for i from 0 <= i < n:
ip = <int>indx.getVal(i,0)
su = b.getVal(ip,0)
b.setVal(ip,0, b.getVal(i,0))
if (ii):
for j from ii <= j < (i-1):
su -= a.getVal(i,j) * b.getVal(j,0)
elif (su):
ii = i
b.setVal(i, 0, su)
for i from n > i >= 0:
su = b.getVal(i,0)
for j from (i+1) <= j < n:
su -= a.getVal(i,j) * b.getVal(j,0)
b.setVal(i, 0, su/a.getVal(i,i))
return b
@cython.boundscheck(False)
@cython.wraparound(False)
cdef MatrixList ludcmp(Matrix a):
#Given a matrix a_{nxn}, this routine replaces it by the LU decomposition of a row-wise permutation of itself.
cdef MatrixList LU = MatrixList()
cdef int n = a.rows
cdef int i, j, k, imax
cdef double big, dum, su, temp
cdef Matrix vv = Matrix(n,1)
cdef Matrix indx = Matrix(n,1) #an output vector that records the row permutation effected by the partial pivoting
cdef Matrix d = Matrix(1,1, ones= True) #an output as +1 or -1 depending on whether the number of row interchanges was even or odd, respectively
cdef double TINY = 1.1e-16
for i from 0 <= i < n:
big = 0.0
for j from 0 <= j < n:
temp=fabs(a.getVal(i,j))
if (temp > big):
big=temp
if (big ==0.0):
raise Exception("ERROR! ludcmp: Singular matrix\n")
vv.setVal(i,0,1.0/big)
for j from 0 <= j < n:
for i from 0 <= i < j:
su = a.getVal(i,j)
for k from 0 <= k < i:
su -= a.getVal(i,k)*a.getVal(k,j)
a.setVal(i,j,su)
big=0.0
for i from j<= i< n:
su = a.getVal(i,j)
for k from 0 <= k < j:
su -= a.getVal(i,k)*a.getVal(k,j)
a.setVal(i, j, su)
dum=vv.getVal(i,0)*fabs(su )
if (dum >= big):
big=dum
imax=i
if (j != imax):
for k from 0 <= k < n:
dum = a.getVal(imax,k)
a.setVal(imax, k, a.getVal(j,k))
a.setVal(j,k, dum)
d.setVal(0, 0, -d.getVal(0,0))
vv.setVal(imax, 0, vv.getVal(j, 0))
indx.setVal(j, 0, imax)
if (a.getVal(j,j) == 0.0):
a.setVal(j,j, TINY)
if (j != (n-1)):
dum=1.0/a.getVal(j,j)
for i from (j+1)<= i <n:
a.setVal(i,j, a.getVal(i,j)*dum)
LU.append(<Matrix>a)
LU.append(<Matrix>indx)
LU.append(<Matrix>d)
return LU
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 lubksb(Matrix a, Matrix indx, Matrix b)
cdef MatrixList ludcmp(Matrix a)
任何帮助查找错误的帮助将不胜感激。
示例：
import numpy
from matrix import Matrix
from numpy.linalg import inv
import timeit
import numpy as np
r=numpy.random.random((100, 100))
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])
start = timeit.default_timer()
x=d.Inv()
stop = timeit.default_timer()
print "LU decomposition:", stop - start
我发现我在 lubksb 函数中犯了一些非常小的错误，通过修复这些错误我得到了正确的答案。这是固定代码：
@cython.boundscheck(False)
@cython.wraparound(False)
cdef Matrix lubksb(Matrix a, Matrix indx, Matrix b):
cdef int n = a.rows
cdef int i, ip, j
cdef int ii = 0
cdef double su
for i from 0 <= i < n:
ip = <int>indx.getVal(i,0)
su = b.getVal(ip,0)
b.setVal(ip,0, b.getVal(i,0))
if (ii>=0):
for j from ii <= j <= (i-1):
su -= a.getVal(i,j) * b.getVal(j,0)
elif (su):
ii = i
b.setVal(i, 0, su)
for i from n >= i >= 0:
su = b.getVal(i,0)
for j from (i+1) <= j < n:
su -= a.getVal(i,j) * b.getVal(j,0)
if (a.getVal(i,i)==0.0):
a.setVal(i,i, 1.1e-16)
b.setVal(i, 0, su/a.getVal(i,i))
return b