基于双边jacobi的奇异值分解Matlab代码

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

function [U, S, V] = bidiag_svd(A)
% Bidiagonalization using Householder reflections
[m, n] = size(A);
B = zeros(m, n);
U = eye(m);
V = eye(n);
for j = 1:n
    [v, beta] = house(A(j:m, j));
    A(j:m, j:n) = A(j:m, j:n) - beta * v * (v' * A(j:m, j:n));
    B(j, j:n) = A(j, j:n);
    if j < n
        [v, beta] = house(A(j:j+1, j+1:n)' );
        A(j:j+1, j+1:n) = A(j:j+1, j+1:n) - beta * v * (v' * A(j:j+1, j+1:n));
        B(j+1:m, j+1) = A(j+1:m, j+1);
    end
    U(j:m, :) = U(j:m, :) - beta * U(j:m, :) * v * v';
end
% SVD of bidiagonal matrix B
[U, S, V] = svd(B);
% Compute left singular vectors of original matrix A
U = U * U';
U = U(:, 1:n);
% Compute right singular vectors of original matrix A
V = V * V';
V = V(:, 1:n);
end

function [v, beta] = house(x)
% Compute Householder reflection vector v and scaling factor beta.
% H = I - beta * v * v'
n = length(x);
sigma = x(2:n)' * x(2:n);
v = [1; x(2:n)];
if sigma == 0 && x(1) >= 0
    beta = 0;
elseif sigma == 0 && x(1) < 0
    beta = -2;
else
    mu = sqrt(x(1)^2 + sigma);
    if x(1) <= 0
        v(1) = x(1) - mu;
    else
        v(1) = -sigma / (x(1) + mu);
    end
    beta = 2 * v(1)^2 / (sigma + v(1)^2);
    v = v / v(1);
end
end

该代码实现了基于双边Jacobi的奇异值分解算法，通过Householder反射将矩阵A转化为上下三角矩阵B，然后对B进行奇异值分解，最后计算出原始矩阵A的左右奇异向量。