线性代数学习四，矩阵的运算和求逆
矩阵的运算是通过Lua的metatable来完成的，用起来和C++的操作符重载一样方便，而矩阵的求逆就是通过求伴随矩阵再除以其行列式的值而得。即
在实现过程中遇到了一些麻烦，有一个是我的对于伴随矩阵的理解有问题，原来伴随矩阵不仅仅要每个元素由其代数余子式代替，还得转置一下。
-
#!/usr/bin/env lua
-
-
require 'determinant'
-
-
function Matrix(data,row,column)
-
local size=table.maxn(data)
-
assert(size==row*column)
-
-
local matrix={}
-
matrix.data=data
-
matrix.row=row
-
matrix.column=column
-
-
matrix.at=function(this,row,column)
-
return this.data[(row-1)*this.column+column]
-
end
-
-
matrix.set=function(this,row,column,num)
-
this.data[(row-1)*this.column+column]=num
-
end
-
-
matrix.swaprow=function(this,r1,r2)
-
local temp
-
for i=1,this.column do
-
temp=this:at(r1,i)
-
this:set(r1,i,this:at(r2,i))
-
this:set(r2,i,temp)
-
end
-
end
-
-
matrix.simplifycolumn=function(this,column,layer)
-
if this:at(layer,column) == 0 then
-
for i=layer+1,this.row do
-
if this:at(i,column) ~= 0 then
-
this:swaprow(i,layer)
-
break
-
end
-
end
-
end
-
-
if this:at(layer,column) == 0 then
-
return layer
-
end
-
-
for i=1,this.row do
-
if i~= layer then
-
local tomultiply=-this:at(i,column)/this:at(layer,column)
-
for j=1,this.column do
-
this:set(i,j,this:at(layer,j)*tomultiply+this:at(i,j))
-
end
-
end
-
end
-
-
return layer+1
-
end
-
-
matrix.swapcolumn=function(this,c1,c2)
-
local temp
-
for i=1,this.row do
-
temp=this:at(i,c1)
-
this:set(i,c1,this:at(i,c2))
-
this:set(i,c2,temp)
-
end
-
end
-
-
matrix.simplifyrow=function(this,row,layer)
-
if this:at(row,layer) == 0 then
-
for i=layer+1,this.column do
-
if this:at(row,i) ~= 0 then
-
this:swapcolumn(i,layer)
-
break
-
end
-
end
-
end
-
-
if this:at(row,layer) == 0 then
-
return layer
-
end
-
-
for i=1,this.column do
-
if i ~= layer then
-
local tomultiply=-this:at(row,i)/this:at(row,layer)
-
for j=1,this.row do
-
this:set(j,i,this:at(j,layer)*tomultiply+this:at(j,i))
-
end
-
end
-
end
-
-
return layer+1
-
end
-
-
-- Make the matrix to reduced form
-
matrix.reduce=function(this)
-
local layer=1
-
for i=1,this.column do
-
layer=this:simplifycolumn(i,layer)
-
if layer>this.row then
-
break
-
end
-
end
-
-
-
for i=1,this.row do
-
local tomultiply=nil
-
for j=1,this.column do
-
if tomultiply then
-
this:set(i,j,this:at(i,j)*tomultiply)
-
elseif this:at(i,j)~= 0 then
-
tomultiply=1/this:at(i,j)
-
this:set(i,j,1)
-
end
-
end
-
end
-
end
-
-
-- Make the matrix to canonical form
-
matrix.canonicalize=function(this)
-
this:reduce()
-
local layer=1
-
for i=1,this.row do
-
layer=this:simplifyrow(i,layer)
-
if layer>this.column then
-
break
-
end
-
end
-
end
-
-
-- Get the rank of the matrix
-
-- Matrix should be reduced before.
-
matrix.rank=function(this)
-
local rank=0
-
for i=1,this.row do
-
for j=1,this.column do
-
if this:at(i,j) ~= 0 then
-
rank=rank+1
-
break
-
end
-
end
-
end
-
-
return rank
-
end
-
-
matrix.print=function(this)
-
for i=1,this.row do
-
for j=1,this.column do
-
io.write(this:at(i,j),'\t')
-
end
-
io.write('\n')
-
end
-
end
-
-
-- Get the transposed matrix
-
matrix.transpose=function(this)
-
local data={}
-
-
for i=1,this.column do
-
for j=1,this.row do
-
table.insert(data,this:at(j,i))
-
end
-
end
-
-
return Matrix(data,this.column,this.row)
-
end
-
-
matrix.cofactor=function(this,row,column)
-
local data={}
-
-
for i=1,this.row do
-
for j=1,this.column do
-
if i~=row and j~=column then
-
table.insert(data,this:at(i,j))
-
end
-
end
-
end
-
-
if (row+column)%2==0 then
-
return Determinant(data):eval()
-
else
-
return -Determinant(data):eval()
-
end
-
end
-
-
-- Get the adjoint matrix of this matrix
-
matrix.adjoint=function(this)
-
local data={}
-
for i=1,this.row do
-
for j=1,this.column do
-
table.insert(data,this:cofactor(i,j))
-
end
-
end
-
-
return Matrix(data,this.row,this.column):transpose()
-
end
-
-
-- Get the inverse matrix of this matrix
-
-- If this matrix is not invertible, nil is returned
-
matrix.inverse=function(this)
-
-- Only array can be invertible
-
if this.row ~= this.column then
-
return nil
-
end
-
-
local det=Determinant(this.data):eval();
-
if det==0 then
-
return nil
-
end
-
-
return (1/det)*this:adjoint()
-
end
-
-
if MatrixCalc then
-
setmetatable(matrix,MatrixCalc)
-
end
-
-
return matrix
-
end
-
-
-- The matrix calculation metatable
-
MatrixCalc={}
-
-
MatrixCalc.__add=function(a,b)
-
assert(a.row==b.row and a.column==b.column)
-
local data={}
-
-
for i=1,a.row*a.column do
-
table.insert(data,a.data[i]+b.data[i])
-
end
-
-
return Matrix(data,a.row,a.column)
-
end
-
-
MatrixCalc.__sub=function(a,b)
-
assert(a.row==b.row and a.column==b.column)
-
local data={}
-
-
for i=1,a.row*a.column do
-
table.insert(data,a.data[i]-b.data[i])
-
end
-
-
return Matrix(data,a.row,a.column)
-
end
-
-
MatrixCalc.__unm=function(a)
-
local data={}
-
-
for i=1,a.row*a.column do
-
table.insert(data,-a.data[i])
-
end
-
-
return Matrix(data,a.row,a.column)
-
end
-
-
MatrixCalc.__mul=function(a,b)
-
if tonumber(a) then
-
-- a is a number
-
local data={}
-
for i=1,b.row*b.column do
-
table.insert(data,a*b.data[i])
-
end
-
return Matrix(data,b.row,b.column)
-
else
-
-- a is a matrix
-
assert(a.column==b.row)
-
-
local data={}
-
for i=1,a.row do
-
for j=1,b.column do
-
local value=0
-
for m=1,a.column do
-
value=value+a:at(i,m)*b:at(m,j)
-
end
-
table.insert(data,value)
-
end
-
end
-
-
return Matrix(data,a.row,b.column)
-
end
-
end
-
-
MatrixCalc.__eq=function(a,b)
-
if a.row~=b.row or a.row~=b.column then
-
return false
-
end
-
-
for i=1,a.row*a.column do
-
if a.data[i]~=b.data[i] then
-
return false
-
end
-
end
-
-
return true
-
end
评论 (0)