33 Julia中的矩阵计算功能

33.1 基本矩阵运算

Julia是一个向Matlab软件一样强调科学计算, 尤其是线性代数运算的编程语言, 所以其矩阵功能比R还要强大。

Julia的基本矩阵构造、加减法、乘法、求逆、求解线性方程组等运算, 见第@(jintro)章。

size(A)返回包含行数和列数的二元组(tuple)。

eye(n)生成n阶的元素类型为Float64的单位阵。 ones(n, m)生成\(n \times m\)的元素类型为Float64的都是1的矩阵, zeros(n,m)生成\(n \times m\)的元素类型为Float64的都是0的矩阵。

A'表示A的共轭转置。 A + BA - B表示矩阵加减, A * B表示矩阵乘法。 inv(A)返回逆矩阵。 A \ b求解\(A x = b\)中的\(x\)

dot(x, y)计算两个向量的内积, 与sum(x .* y)相同。 norm(x)求向量的欧式模, 等于sqrt(sum(x .^ 2))

factorize(A)A计算适合的分解, 自动判断A的特殊结构。

det(A)求行列式。 trace(A)计算迹。

eigvals(A)求特征值。 eig(A)返回由特征值、特征向量组成的二元组, 特征值储存在一个向量中, 特征向量储存在矩阵的各列中。

eig(A, B)求解广义特征值问题\(A v = \lambda B v\)

rank(A)求秩。

Julia的线性代数模块提供了丰富的矩阵类型和矩阵计算功能, 详见Base.LinAlg模块的文档。

33.2 用Givens变换作QR分解示例

using Distributions # 随机数函数
srand(101) # 指定起始种子
A=round.(rand(4,4), 2)

##  4×4 Array{Float64,2}:
##   0.09  1.0   0.97  0.42
##   0.83  0.86  0.62  0.8 
##   0.63  0.8   0.75  0.23
##   0.16  0.95  0.73  0.07

## Givens变换矩阵,在i1,i2两行和两列作用,将(a, b)长度集中到a所在元素
function givensmat(n, a, b, i1, i2)
    submat = [a b; -b a]/sqrt(a*a+b*b)
    resmat = eye(n)
    resmat[[i1,i2], [i1,i2]] = submat
    resmat
end

n=4
Q1 = givensmat(n, A[1,1], A[2,1], 1, 2)

##     4×4 Array{Float64,2}:
##       0.107802  0.994172  0.0  0.0
##      -0.994172  0.107802  0.0  0.0
##       0.0       0.0       1.0  0.0
##       0.0       0.0       0.0  1.0

A12 = Q1 * A

##     4×4 Array{Float64,2}:
##      0.834865      0.96279    0.720955   0.840615
##      1.38778e-17  -0.901463  -0.89751   -0.331311
##      0.63          0.8        0.75       0.23    
##      0.16          0.95       0.73       0.07    

Q2 = givensmat(n, A12[1,1], A12[3,1], 1, 3)

##     4×4 Array{Float64,2}:
##       0.798229  0.0  0.602354  0.0
##       0.0       1.0  0.0       0.0
##      -0.602354  0.0  0.798229  0.0
##       0.0       0.0  0.0       1.0

A13 = Q2 * A12

##     4×4 Array{Float64,2}:
##      1.0459        1.25041     1.02725    0.809545
##      1.38778e-17  -0.901463   -0.89751   -0.331311
##      0.0           0.0586429   0.164402  -0.322755
##      0.16          0.95        0.73       0.07    

Q3 = givensmat(n, A13[1,1], A[4,1], 1, 4)

##     4×4 Array{Float64,2}:
##       0.9885   0.0  0.0  0.15122
##       0.0      1.0  0.0  0.0    
##       0.0      0.0  1.0  0.0    
##      -0.15122  0.0  0.0  0.9885 

A14 = Q3 * A13

##     4×4 Array{Float64,2}:
##       1.05806       1.37969     1.12583    0.81082  
##       1.38778e-17  -0.901463   -0.89751   -0.331311 
##       0.0           0.0586429   0.164402  -0.322755 
##      -2.77556e-17   0.749989    0.566264  -0.0532239

Q4 = givensmat(n, A14[2,2], A14[3,2], 2, 3)

##     4×4 Array{Float64,2}:
##      1.0   0.0         0.0        0.0
##      0.0  -0.997891    0.0649159  0.0
##      0.0  -0.0649159  -0.997891   0.0
##      0.0   0.0         0.0        1.0

A23 = Q4 * A14

##     4×4 Array{Float64,2}:
##       1.05806       1.37969       1.12583    0.81082  
##      -1.38485e-17   0.903368      0.906289   0.30966  
##      -9.00889e-19  -6.93889e-18  -0.105793   0.343581 
##      -2.77556e-17   0.749989      0.566264  -0.0532239

Q5 = givensmat(n, A23[2,2], A23[4,2], 2, 4);
A24 = Q5 * A23

##     4×4 Array{Float64,2}:
##       1.05806       1.37969       1.12583    0.81082 
##      -2.83844e-17   1.17412       1.05901    0.204255
##      -9.00889e-19  -6.93889e-18  -0.105793   0.343581
##      -1.25092e-17   0.0          -0.143223  -0.238751

Q6 = givensmat(n, A24[3,3], A24[4,3], 3, 4);
A34 = Q6 * A24

##     4×4 Array{Float64,2}:
##       1.05806       1.37969      1.12583    0.81082 
##      -2.83844e-17   1.17412      1.05901    0.204255
##       1.05971e-17   4.1227e-18   0.178059  -0.012095
##       6.70761e-18  -5.58136e-18  0.0        0.418215

到A34已经变成了上三角矩阵。各个Givens旋转合并起来也是一个正交阵:

上面相当于\(Q A = R\),所以\(A = Q' R\)。 Julia的qr()函数可以计算QR分解:

qr()函数给出的上三角矩阵不满足对角线元素为正的要求。需要时可以将负的主对角元素所在行反号并将对应的Q的列反号。如: