34 Julia中的矩阵计算功能(*)
34.1 基本矩阵运算
Julia是一个像Matlab软件一样强调科学计算， 尤其是线性代数运算的编程语言， 所以其矩阵功能比R还要强大。
Julia的基本矩阵构造、加减法、乘法、求逆、求解线性方程组等运算， 见附录B。
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 + B和A - 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模块的文档。
34.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 = Q6*Q5*Q4*Q3*Q2*Q1
## 4×4 Array{Float64,2}:
## 0.085061 0.784451 0.595427 0.15122
## 0.751748 -0.189333 -0.0183152 0.631421
## 0.43877 -0.351868 0.556258 -0.611757
## 0.484892 0.474318 -0.579403 -0.451877
Q'*Q
## 4×4 Array{Float64,2}:
## 1.0 -5.55112e-17 -5.55112e-17 0.0
## -5.55112e-17 1.0 -2.22045e-16 0.0
## -5.55112e-17 -2.22045e-16 1.0 0.0
## 0.0 0.0 0.0 1.0
B = Q * A
## 4×4 Array{Float64,2}:
## 1.05806 1.37969 1.12583 0.81082
## 0.0 1.17412 1.05901 0.204255
## 4.16334e-17 0.0 0.178059 -0.012095
## -4.16334e-17 5.55112e-17 1.11022e-16 0.418215
maximum(B - A34)
## 2.220446049250313e-16上面相当于\(Q A = R\)，所以\(A = Q' R\)。
Julia的qr()函数可以计算QR分解：
Q1, R1 = qr(A);
R1
## 4×4 Array{Float64,2}:
## -1.05806 -1.37969 -1.12583 -0.81082
## 0.0 1.17412 1.05901 0.204255
## 0.0 0.0 -0.178059 0.012095
## 0.0 0.0 0.0 -0.418215qr()函数给出的上三角矩阵不满足对角线元素为正的要求。需要时可以将负的主对角元素所在行反号并将对应的Q的列反号。如：
R2 = R1; Q2 = Q1
R2[1,:] = -R2[1,:]; Q2[:,1] = -Q2[:,1]
R2[3,:] = -R2[3,:]; Q2[:,3] = -Q2[:,3]
R2[4,:] = -R2[4,:]; Q2[:,4] = -Q2[:,4]
R2
## 4×4 Array{Float64,2}:
## 1.05806 1.37969 1.12583 0.81082
## 0.0 1.17412 1.05901 0.204255
## -0.0 -0.0 0.178059 -0.012095
## -0.0 -0.0 -0.0 0.418215
Q2
## 4×4 Array{Float64,2}:
## 0.085061 0.751748 0.43877 0.484892
## 0.784451 -0.189333 -0.351868 0.474318
## 0.595427 -0.0183152 0.556258 -0.579403
## 0.15122 0.631421 -0.611757 -0.451877
Q'
## 4×4 Array{Float64,2}:
## 0.085061 0.751748 0.43877 0.484892
## 0.784451 -0.189333 -0.351868 0.474318
## 0.595427 -0.0183152 0.556258 -0.579403
## 0.15122 0.631421 -0.611757 -0.451877