基于双边jacobi的svd算法Matlab代码和数值模拟

双边Jacobi的SVD算法是一种经典的求解SVD分解的方法，其主要思想是通过迭代的方式将矩阵分解成奇异值、左奇异向量和右奇异向量的乘积形式。下面给出Matlab代码和数值模拟。

Matlab代码：

function [U,S,V] = bidiag_jacobi_svd(A)
% Bidiagonalization with Jacobi rotations
% Input: A (m x n) matrix
% Output: U, S, V such that A = U*S*V'
[m,n] = size(A);
U = eye(m);
V = eye(n);
B = A;
for k = 1:n
    % Compute Householder reflection to zero subdiagonal entries
    [v,b] = house(B(k:m,k));
    V_k = eye(n);
    V_k(k:n,k:n) = eye(n-k+1) - b*v*v';
    B = V_k*B;
    % Compute Householder reflection to zero superdiagonal entries
    if k < n
        [v,b] = house(B(k,k+1:n)');
        U_k = eye(m);
        U_k(k:m,k:m) = eye(m-k+1) - b*v*v';
        B = B*U_k';
        V(:,k+1:n) = V(:,k+1:n)*U_k;
    end
end
% Apply QR algorithm to diagonal and superdiagonal entries
for k = n:-1:1
    while norm(B(k,k+1:n)) > eps*norm(B)
        % Compute Wilkinson shift
        mu = wilkinson_shift(B(k-1:k,k-1:k));
        % Apply shifted QR iteration
        [Q,R] = qr(B(1:k,1:k) - mu*eye(k));
        B(1:k,1:k) = R*Q + mu*eye(k);
        U(:,1:k) = U(:,1:k)*Q;
        V(:,1:k) = V(:,1:k)*Q;
    end
end
S = diag(B);
end

function [v,b] = house(x)
% Compute Householder reflection to zero subdiagonal entry
n = length(x);
sigma = x(2:n)'*x(2:n);
v = [1; x(2:n)];
if sigma == 0
    b = 0;
else
    mu = sqrt(x(1)^2+sigma);
    if x(1) <= 0
        v(1) = x(1) - mu;
    else
        v(1) = -sigma/(x(1)+mu);
    end
    b = 2*v(1)^2/(sigma+v(1)^2);
    v = v/v(1);
end
end

function mu = wilkinson_shift(A)
% Compute Wilkinson shift
n = size(A,1);
d = (A(n-1,n-1)-A(n,n))/2;
mu = A(n,n) - sign(d)*A(n-1,n)^2/(abs(d)+sqrt(d^2+A(n-1,n)^2));
end

数值模拟：

为了验证双边Jacobi的SVD算法的正确性，我们可以使用随机生成的矩阵进行数值模拟。下面是一个简单的示例代码：

n = 100;
m = 50;
A = randn(m,n);
[U,S,V] = bidiag_jacobi_svd(A);
subplot(2,1,1);
semilogy(diag(S),'o');
title('Singular values');
subplot(2,2,3);
imagesc(A);
title('Original matrix');
subplot(2,2,4);
imagesc(U*S*V');
title('Reconstructed matrix');

这段代码会生成一个 $50 \times 100$ 的随机矩阵，并对其进行SVD分解和重构。我们可以使用 `semilogy` 函数绘制奇异值的对数图像，以及使用 `imagesc` 函数绘制原始矩阵和重构矩阵的色块图像。如果算法正确，则奇异值应该是单调非增的，重构矩阵应该与原始矩阵非常接近。