这学期刚学完数值代数,在这里总结一下学过的算法,附MATLAB代码。
QR分解有两种形式:完全QR分解(reduced QR factorization)、约化QR分解(full QR factorization)。
完全QR分解:,其中
,
为酉矩阵,
为上三角矩阵。
约化QR分解:,其中
,
具有正交列,
为上三角矩阵。
实际上,两者是等价的。因为比
多的几列可以通过扩充标准正交基得到,
比
多的几行全为0。我认为完全QR分解最大的优点就是
是酉矩阵,很好处理。
以下是三中QR分解的经典算法,其中豪斯霍尔德算法工作量最少且向后稳定。
1.经典格拉姆-施密特正交化(classical Gram-Schmidt iteration)
基本思想:
找到单位正交序列满足
这等价于
利用这种思想对矩阵进行约化QR分解。
代码实现
function [Q,R]=classical_gramschmidt(A)
[~,n] = size(A);
V = zeros(size(A));
Q = zeros(size(A));
R = zeros(n);
for j = 1:n
V(:,j) = A(:,j);
for i = 1:j-1
R(i,j) = Q(:,i)'*A(:,j);
V(:,j) = V(:,j)-R(i,j)*Q(:,i);
end
R(j,j) = norm(V(:,j),2);
Q(:,j) = V(:,j)/R(j,j);
end
end
2.修正的格拉姆-施密特正交化(modified Gram-Schmidt iteration)
基本思想
设为投影到正交于
的正交投影算子,这样做可以保证
,
由正交投影算子的相关知识可以得到,,其中
代码实现
function [Q,R]=modified_gramschmidt(A)
[~,n] = size(A);
Q = zeros(size(A));
R = zeros(n);
V=A;
for i = 1:n
R(i,i) = norm(V(:,i),2);
Q(:,i) = V(:,i)/R(i,i);
for j = i+1:n
R(i,j) = Q(:,i)'*V(:,j);
V(:,j) = V(:,j)-R(i,j)*Q(:,i);
end
end
end
3.豪斯霍尔德三角形化(Householder)
function [Q,A] = householder(A)
[m,n] = size(A);
Q = eye(m);
for k = 1:n
x = A(k:m,k);
v = sign(x(1))*norm(x,2)*eye(m-k+1,1)+x;
v = v/norm(v,2);
F = eye(m-k+1)-2*v*v'/(v'*v);
Q = [eye(k-1) zeros(k-1,m-k+1);zeros(m-k+1,k-1) F]*Q;
A(k:m,k:n) = A(k:m,k:n)-2*v*(v'*A(k:m,k:n));
end
Q = Q';
end