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

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

function [U, S, V] = bidiagonal_jacobi(A)
% Bidiagonal Jacobi SVD algorithm
% A: m x n matrix
% U: m x m orthogonal matrix
% S: m x n diagonal matrix with singular values
% V: n x n orthogonal matrix

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

for i = 1:min(m, n)
    % Bidiagonalization
    [U, S, V] = bidiag_jacobi_step(U, S, V, i);
end

% Truncate small values in S
tol = max(m, n) * eps(norm(S, inf));
S(abs(S) < tol) = 0;

% Sort singular values and columns of U and V accordingly
[S, idx] = sort(diag(S), 'descend');
U = U(:, idx);
V = V(:, idx);
end

function [U, S, V] = bidiag_jacobi_step(U, S, V, k)
% Perform one step of bidiagonal Jacobi algorithm
% on S(k:m, k:n)

[m, n] = size(S);
sigma = S(k, k);
beta = norm(S(k:end, k));

if beta == 0
    % S(k:end, k) is already zero
    return
end

gamma = sqrt(sigma^2 + beta^2);
if sigma < 0
    gamma = -gamma;
end

c = sigma / gamma;
s = beta / gamma;

% Apply Givens rotations to S, U, and V
G1 = [c, -s; s, c];
S(k:k+1, k:n) = G1 * S(k:k+1, k:n);
if k < m
    G2 = [c, -s; s, c];
    S(k:m, k:k+1) = S(k:m, k:k+1) * G2';
    U(:, k:k+1) = U(:, k:k+1) * G2';
end
if k < n
    G3 = [c, -s; s, c];
    S(k:k+1, k+1:n) = G3' * S(k:k+1, k+1:n);
    V(:, k:k+1) = V(:, k:k+1) * G3;
end
end

示例：

% Generate a random m x n matrix
m = 5;
n = 4;
A = rand(m, n);

% Compute SVD using built-in function for comparison
[U0, S0, V0] = svd(A);

% Compute SVD using bidiagonal Jacobi algorithm
[U, S, V] = bidiagonal_jacobi(A);

% Check if the result is close to the built-in function
fprintf('Maximum absolute difference in U: %g\n', max(abs(U - U0), [], 'all'));
fprintf('Maximum absolute difference in S: %g\n', max(abs(S - S0), [], 'all'));
fprintf('Maximum absolute difference in V: %g\n', max(abs(V - V0), [], 'all'));