基于双边Jacobi的奇异值分解的Matlab代码和示例

下面是基于双边Jacobi的奇异值分解的Matlab代码和示例：

function [U, S, V] = bidiagjacobi(A)
% Bidiagonalization using Jacobi rotation
% A: input matrix
% U: left orthogonal matrix
% S: bidiagonal matrix with singular values on diagonal
% V: right orthogonal matrix

[m, n] = size(A);
U = eye(m);
V = eye(n);
S = A;

for k = 1:n
    % Perform Jacobi rotations to zero out entries below diagonal
    for i = k+1:m
        theta = atan2(S(i,k), S(k,k));
        c = cos(theta);
        s = sin(theta);
        G = [c -s; s c];
        S([k,i],k:n) = G' * S([k,i],k:n);
        S(k:n,[k,i]) = S(k:n,[k,i]) * G;
        U([k,i],:) = G' * U([k,i],:);
    end
    
    % Perform Jacobi rotations to zero out entries to the right of diagonal
    if k < n
        for j = k+1:n
            theta = atan2(S(k,j), S(k,k));
            c = cos(theta);
            s = sin(theta);
            G = [c -s; s c];
            S(k:m,[k,j]) = S(k:m,[k,j]) * G;
            S(:,[k,j]) = S(:,[k,j]) * G';
            V([k,j],:) = G' * V([k,j],:);
        end
    end
end

示例：

% Generate a random matrix
A = rand(5);

% Compute SVD using bidiagonalization with Jacobi rotation
[U, S, V] = bidiagjacobi(A);

% Compare with SVD computed using built-in Matlab function
[Uref, Sref, Vref] = svd(A);
disp(norm(U - Uref))
disp(norm(S - Sref))
disp(norm(V - Vref))

运行示例代码会输出三个值，分别是 $U$ 和 $U_{ref}$ 之间的 Frobenius 范数，$S$ 和 $S_{ref}$ 之间的 Frobenius 范数，以及 $V$ 和 $V_{ref}$ 之间的 Frobenius 范数。这些值都应该非常接近 0，证明我们的实现和 Matlab 内置的 SVD 函数得到了相同的结果。