编程用矩阵运算解线性方程组

为什么80%的码农都做不了架构师？>>>   

矩阵的乘法：

2	3		3
		×		=	?
3	5		4 

ruby标准库中有Matrix，定义矩阵是Matrix[]形式：

2.4.0 :001 > require 'matrix'
 => true 
2.4.0 :002 > Matrix[[2, 3], [3, 5]] * Matrix[[3],[4]]
 => Matrix[[18], [29]] 

Julia：可以直接以数组形式来写矩阵，多行用分号分隔，每列用空格分隔。

julia> [2 3; 3 5] * [3; 4] 
2-element Array{Int64,1}:
 18
 29

现在：记2元线性方程组：

2x + 3y = 18
3x + 5y = 29

可以抽象成矩阵乘法为：

2	3		x		18
		×		=		
3	5 		y		29

即 Matrix[[2,3],[3,5]] × Matrix[[x],[y]] = Matrix[[18],[29]] Ruby的标准库中有直接求逆矩阵：A.inv。 定义 A × B = C，根据逆矩阵变换，可以推知：B = A.inv × C B为什么是A的逆矩阵 × C呢？这里可以简单证明：

A × B = C
A逆 × A × B = A逆 × C
A逆 × A = I (单位矩阵）
I × B = A逆 × C

I可以直接消除(单位矩阵乘以任何矩阵都不影响矩阵）(I类似于数值运算中的1，即 1 × n = n 和 n × 1 = n）。 故：B = A逆 × C julia中有一个大写的字母I直接表示任何长度的单位矩阵（意味着他能和任意矩阵相乘） I定义：

An object of type UniformScaling, representing an identity matrix of any size.

julia> I
UniformScaling{Int64}
1*I

julia> I * [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> [1 2; 3 4] * I
2×2 Array{Int64,2}:
 1  2
 3  4

通过Julia的语言验证，你应该可以确信，上面的I可以直接消除了吧。

现在求[x,y]的值：

2.4.0 :003 > Matrix[[2,3],[3,5]].inv * Matrix[[18],[29]]
 => Matrix[[(3/1)], [(4/1)]] 

结果为分数形式：3/1和4/1即 3,4。

Julia： Julia中，定义 A × B = C 则 A \ C = B （ \ 为矩阵左除运算符，A \ B = inv(A) * B）

julia> [2 3; 3 5] \ [18; 29]
2-element Array{Float64,1}:
 3.0
 4.0

当然了，也可以用逆矩阵来解：

julia> inv([2 3; 3 5]) * [18; 29]
2-element Array{Float64,1}:
 3.0
 4.0

相比较Ruby，我更喜欢Julia的矩阵除法操作，定义和运算都很简洁。 Julia还有一个矩阵的右除运算符（/，就是标准的除法，A / B = A * inv(B) )，这个右除即等同于Ruby的矩阵 / 方法。

由上面可知道，如果一个语言内部库实现了逆矩阵以及矩阵乘法，就可以用来解线性方程组了。（当然有些逆矩阵是不存在，排除掉这种情况。）

Ruby中虽然也定义了一个矩阵的除法，不过其除法的定义是： A / B = A × B.inverse

而矩阵的乘法是不满足交换律的：即 A × B != B × A 所以，无法通过Ruby的矩阵除法来求解上面格式的线性方程组，只能通过乘法来求： 即 A.inv × C = B

当然了，如果要用这个除法来解（Ruby的矩阵除方法 / 和 Julia的矩阵右除函数 / ） ，我们只能将线性方程做一次变换，即写成：B x A = C，这样 B = C / A 了：

然后我们就可以解了：

2.4.0 :004 > Matrix[[18,29]] / Matrix[[2,3],[3,5]]
 => Matrix[[(3/1), (4/1)]] 

Julia：

julia> [18 29] / [2 3; 3 5]
1×2 Array{Float64,2}:
3.0  4.0

Wolframalpha： 对于矩阵右除，只要会正常改写线性方程组为矩阵乘法形式，不会很复杂。

Haskell原生没有矩阵库，不过可以通过cabal来安装matrix： https://hackage.haskell.org/package/matrix-0.3.5.0/docs/Data-Matrix.html

cabal update
cabal install matrix
ghci -package matrix

进入ghci后，加载 matrix：位于Data.Matrix

Prelude> :m Data.Matrix
Prelude Data.Matrix> fromLists [[2,3],[3,5]] * fromLists [[3],[4]]
( 18 )
( 29 )

写一个函数来实现类似于julia的矩阵除法 （这里用双反斜杠来表示）

Prelude Data.Matrix> a \\ b = (flip (*) b) <$> (inverse a)
Prelude Data.Matrix> fromLists [[2,3],[3,5]] \\ fromLists [[18],[29]]
Right ( 3.0 )
( 4.0 )

可以看到，正常求出了x,y为3,4 定义的函数为Matrix类：

Prelude Data.Matrix> :t (\\)
(\\)
  :: (Fractional a, Eq a) =>
     Matrix a -> Matrix a -> Either String (Matrix a)

也可以定义一个List类上的运算：（这样就更像Julia的调用方式了）

Prelude Data.Matrix> a \\\ b = (flip (*) $ fromLists b) <$> (inverse $ fromLists a)
Prelude Data.Matrix> [[2,3],[3,5]] \\\ [[18],[29]]
Right ( 3.0 )
( 4.0 )

查看类型：

(\\\)
  :: (Fractional a, Eq a) =>
     [[a]] -> [[a]] -> Either String (Matrix a)

实现Haskell的矩阵右除法：B X A = C 得 B = C / A ：

Prelude Data.Matrix> let a / b  = (a*) <$> (inverse b)
Prelude Data.Matrix> fromLists [[18,29]] / fromLists [[2,3],[3,5]]
Right ( 3.0 4.0 )

转载于:https://my.oschina.net/jsk/blog/911115