Function svd 

pub fn svd<const D: usize, const D1: usize>(
    tensor: Tensor<D>,
    sweeps: usize,
) -> (Tensor<D>, Tensor<D1>, Tensor<D>)

Computes the singular value decomposition of a square or rectangular matrix.

This function decomposes the input tensor A into three tensors U, S, Vt such that A = U @ diag(S) @ Vt, where:

• U has orthonormal columns of shape [..., m, k]
• S contains the singular values of A, sorted in descending order, of shape [..., k]
• Vt has orthonormal rows of shape [..., k, n]

with k = min(m, n) (reduced decomposition, matching torch.linalg.svd).

Algorithm

Two stages, mirroring the LAPACK gesvd / dbdsqr structure:

1. Golub-Kahan bidiagonalization using Householder reflections (A = U1 B V1^T).
2. Implicitly shifted bidiagonal QR iteration (LAPACK dbdsqr) to diagonalize B.

Arguments

• tensor - The input tensor of shape [..., m, n].
• sweeps - Upper bound on the number of QR sweeps per singular value.

Returns

A tuple of three tensors (U, S, Vt):

• U: [..., m, k] with orthonormal columns
• S: [..., k] singular values in descending order
• Vt: [..., k, n] with orthonormal rows

Generic Parameters

• D: The number of dimensions of the input tensor.
• D1: Must be set to D - 1 (the rank of the singular value tensor).

Panics

This function will panic if the tensor checks fail:

• The input tensor has less than 2 dimensions (D < 2).
• The input is a quantized tensor with dtype DType::QFloat.
• The generic parameters do not satisfy D - 1 == D1.
• The input tensor requires gradients (SVD has no autodiff support yet; detach first).
• sweeps is zero.
• The QR iteration does not converge within the requested sweep budget.

Performance Note

The computation is dispatched to the backend through LinalgOps::svd, which backends may override with a native or fused implementation (none ship one yet). The default implementation runs the reference pipeline on the host over the tensor data (into_data / from_data), which is deterministic and backend-independent, but the bidiagonalization is O(m n^2) scalar math. It is not competitive with tuned native libraries (e.g. cuSOLVER) for large matrices.

Numerical Behavior

• If the input tensor has dtype F16 or BF16, it is internally upcast to F32 for the computation and cast back to the original dtype before returning, like det and lu.
• Singular values are sorted in descending order.
• Internal norms, rotations, and shifts use scaled formulations to reduce overflow and underflow for extreme finite inputs.

Example

use burn_linalg::svd;
use burn::Tensor;

fn example() {
    let device = Default::default();
    let tensor = Tensor::<2>::from_data([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], &device);

    // Compute the singular value decomposition
    let (u, s, vt) = svd::<2, 1>(tensor, 10);

    // A = U @ diag(S) @ Vt (within tolerance)
    let recon = u.mul(s.unsqueeze_dim(0)).matmul(vt);
    println!("{}", recon);
}