SingularValueDecomposition.f90

Go to the documentation of this file.

module svdmodule
  use kindmodule, only: dp, lgp, i4b
  use constantsmodule, only: done, dprec, dsame, dem6
  use linearalgebrautilsmodule, only: eye, outer_product
 
  implicit none
  private
 
  public :: svd
  public :: bidiagonal_decomposition
  public :: bidiagonal_qr_decomposition
 
contains
 
  !!> @brief Constructs a Householder reflection matrix for a given vector.
  !!
  !<
  pure function house_holder_matrix(x) result(Q)
    ! -- dummy
    real(dp), intent(in) :: x(:)
    real(dp), allocatable, dimension(:, :) :: q
    ! -- locals
    real(dp) :: x_norm, y_norm
    real(dp), allocatable :: v(:)
    real(dp), allocatable :: w(:)
 
    x_norm = norm2(x)
    y_norm = norm2(x(2:))
 
    q = eye(size(x))
 
    if (dabs(y_norm) < dsame) then
      return
    end if
 
    v = x
    if (x(1) > 0) then
      ! Parlett (1971) suggested this modification to avoid cancellation
      v(1) = -(y_norm**2) / (x(1) + x_norm)
    else
      v(1) = v(1) - x_norm
    end if
 
    w = v / norm2(v)
 
    q = q - 2.0_dp * outer_product(w, w)
 
  end function house_holder_matrix
 
  !!> @brief Reduces a matrix to bidiagonal form using Householder transformations.
  !!
  !! This subroutine transforms the input matrix `A` into a bidiagonal matrix by applying a sequence of
  !! Householder reflections from the left and right. The orthogonal transformations are accumulated in
  !! the matrices `P` (left transformations) and `Qt` (right transformations).
  !!
  !! The result is:
  !!     P^T * A * Qt = B
  !! where `B` is bidiagonal, `P` and `Qt` are orthogonal matrices.
  !!
  !! This step is the first phase of the Golub-Reinsch SVD algorithm and prepares the matrix for further
  !! diagonalization using the QR algorithm.
  !!
  !<
  pure subroutine bidiagonal_decomposition(A, P, Qt)
    ! -- dummy
    real(dp), intent(inout), dimension(:, :) :: a
    real(dp), intent(out), dimension(:, :), allocatable :: p, qt
    ! -- locals
    integer(I4B) :: m, n
    integer(I4B) :: i
    real(dp), allocatable, dimension(:, :) :: qi, pi
    real(dp), allocatable, dimension(:) :: h
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    qt = eye(n)
    p = eye(m)
 
    do i = 1, min(m, n)
      ! columns
      h = a(i:m, i)
      pi = house_holder_matrix(h)
      a(i:m, :) = matmul(pi, a(i:m, :)) ! Apply householder transformation from left
      p(:, i:m) = matmul(p(:, i:m), pi)
 
      ! rows
      if (i < n) then
        h = a(i, i + 1:n)
        qi = transpose(house_holder_matrix(h))
        a(:, i + 1:n) = matmul(a(:, i + 1:n), qi) ! Apply householder transformation from right
        qt(i + 1:n, :) = matmul(qi, qt(i + 1:n, :))
      end if
 
    end do
 
  end subroutine bidiagonal_decomposition
 
  !!> @brief Constructs a Givens rotation matrix.
  !!
  !<
  pure function givens_rotation(a, b) result(G)
    ! -- dummy
    real(dp), intent(in) :: a, b
    real(dp), dimension(2, 2) :: g
    ! -- locals
    real(dp) :: c, s, h, d
 
    if (abs(b) < dprec) then
      g = eye(2)
      return
    end if
 
    h = hypot(a, b)
    d = 1.0 / h
    c = abs(a) * d
    s = sign(d, a) * b
 
    g(1, 1) = c
    g(1, 2) = s
    g(2, 1) = -s
    g(2, 2) = c
 
  end function givens_rotation
 
  !!> @brief Computes the Wilkinson shift for the lower-right 2x2 block of a bidiagonal matrix.
  !!
  !! This function calculates the Wilkinson shift, which is used to accelerate convergence in the implicit QR algorithm
  !! for bidiagonal matrices during SVD. The shift is computed from the lower-right 2x2 block of the input matrix `A`
  !! and is chosen to be the eigenvalue of that block which is closer to the bottom-right element.
  !!
  !! The shift is used to improve numerical stability and speed up the annihilation of off-diagonal elements.
  !!
  !<
  pure function compute_shift(A) result(mu)
    ! -- dummy
    real(dp), intent(in), dimension(:, :) :: a
    real(dp) :: mu
    ! -- locals
    integer(I4B) :: m, n, min_mn
    real(dp) t11, t12, t21, t22
    real(dp) dm, fmmin, fm, dn
    real(dp) :: mean, product, mean_product, mu1, mu2
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    min_mn = min(m, n)
 
    dn = a(min_mn, min_mn)
    dm = a(min_mn - 1, min_mn - 1)
    fm = a(min_mn - 1, min_mn)
    if (min_mn > 2) then
      fmmin = a(min_mn - 2, min_mn - 1)
    else
      fmmin = 0.0_dp
    end if
 
    t11 = dm**2 + fmmin**2
    t12 = dm * fm
    t21 = t12
    t22 = fm**2 + dn**2
 
    mean = (t11 + t22) / 2.0_dp
    product = t11 * t22 - t12 * t21
    mean_product = mean**2 - product
    if (abs(mean_product) < dem6) then
      mean_product = 0.0_dp
    end if
    mu1 = mean - sqrt(mean_product)
    mu2 = mean + sqrt(mean_product)
    if (abs(t22 - mu1) < abs(t22 - mu2)) then
      mu = mu1
    else
      mu = mu2
    end if
 
  end function compute_shift
 
  !!> @brief Performs QR iterations on a bidiagonal matrix as part of the SVD algorithm.
  !!
  !<
  pure subroutine bidiagonal_qr_decomposition(A, U, VT)
    ! -- dummy
    real(dp), intent(inout), dimension(:, :) :: a
    real(dp), intent(inout), dimension(:, :) :: u, vt
    ! -- locals
    integer(I4B) :: m, n, i, j
    real(dp), dimension(2, 2) :: g
    real(dp) :: t11, t12, mu
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    DO i = 1, n - 1
      j = i + 1
      if (i == 1) then
        ! For the first iteration use the implicit shift
        mu = compute_shift(a)
        ! Apply the shift to the full tri-diagonal matrix.
        t11 = a(1, 1)**2 - mu
        t12 = a(1, 1) * a(1, 2)
 
        g = givens_rotation(t11, t12)
      else
        g = givens_rotation(a(i - 1, i), a(i - 1, j))
      end if
 
      a(:, i:j) = matmul(a(:, i:j), transpose(g))
      vt(i:j, :) = matmul(g, vt(i:j, :))
 
      if (j > m) cycle
      g = givens_rotation(a(i, i), a(j, i))
      a(i:j, :) = matmul(g, a(i:j, :))
      u(:, i:j) = matmul(u(:, i:j), transpose(g))
    END DO
 
  end subroutine bidiagonal_qr_decomposition
 
  !!> @brief Handles zero (or near-zero) diagonal elements in a bidiagonal matrix during SVD.
  !!
  !! This subroutine detects and processes zero (or near-zero) diagonal elements in the bidiagonal matrix `A`.
  !! If a zero is found on the diagonal, it applies a sequence of Givens rotations to either zero out the
  !! corresponding column (if the zero is at the end of the diagonal) or the corresponding row (otherwise).
  !! The orthogonal transformations are accumulated in the `U` and `VT` matrices as appropriate.
  !! This step is essential for splitting the matrix into smaller blocks and ensuring numerical stability
  !! during the iterative QR step of the SVD algorithm.
  !!
  !<
  pure subroutine handle_zero_diagonal(A, U, VT)
    ! -- dummy
    real(dp), intent(inout), dimension(:, :) :: a
    real(dp), intent(inout), dimension(:, :) :: u, vt
    ! -- locals
    integer(I4B) :: m, n, i
    real(dp), dimension(2, 2) :: g
    integer(I4B) :: zero_index
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    do i = 1, min(m, n)
      if (abs(a(i, i)) < dprec) then
        zero_index = i
        exit
      end if
    end do
 
    if (zero_index == min(n, m)) then
      ! If the zero index is the last element of the diagonal then zero out the column
      do i = zero_index - 1, 1, -1
        g = givens_rotation(a(i, i), a(i, zero_index))
        a(:, [i, zero_index]) = matmul(a(:, [i, zero_index]), transpose(g))
        vt([i, zero_index], :) = matmul(g, vt([i, zero_index], :))
      end do
    else
      ! Else zero out the row
      do i = zero_index + 1, n
        g = givens_rotation(a(i, i), a(zero_index, i))
        a([zero_index, i], :) = matmul(transpose(g), a([zero_index, i], :))
        u(:, [zero_index, i]) = matmul(u(:, [zero_index, i]), g)
      end do
    end if
  end subroutine handle_zero_diagonal
 
  !!> @brief Computes the infinity norm of the superdiagonal elements of a matrix.
  !!
  !<
  pure function superdiagonal_norm(A) result(norm)
    ! -- dummy
    real(dp), intent(in) :: a(:, :)
    real(dp) :: norm
    ! -- locals
    integer(I4B) :: m, n, i
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    norm = 0.0_dp
    do i = 1, min(m, n) - 1
      norm = max(norm, abs(a(i, i + 1)))
    end do
 
  end function superdiagonal_norm
 
  !!> @brief Finds the range of nonzero superdiagonal elements in a bidiagonal matrix.
  !!
  !! This subroutine scans the superdiagonal of the matrix `A` and determines the indices `p` and `q`
  !! such that all superdiagonal elements outside the range `A(p:q, p:q)` are (near-)zero, while
  !! those within the range may be nonzero. This is useful for identifying active submatrices during
  !! the iterative QR step of the SVD algorithm.
  !!
  !! The search starts from the lower-right corner and moves upward to find the last nonzero
  !! superdiagonal (sets `q`), then moves upward again to find the first zero (sets `p`).
  !<
  pure subroutine find_nonzero_superdiagonal(A, p, q)
    ! -- dummy
    real(dp), intent(in), dimension(:, :) :: a
    integer(I4B), intent(out) :: p ! Index of the first nonzero superdiagonal (start of active block)
    integer(I4B), intent(out) :: q ! Index of the last nonzero superdiagonal (end of active block)
    ! -- locals
    integer(I4B) :: m, n, j, min_mn
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    min_mn = min(m, n)
    p = 1
    q = min_mn
 
    do j = min_mn, 2, -1
      if (abs(a(j - 1, j)) > dprec) then
        q = j
        exit
      end if
    end do
 
    do j = q - 1, 2, -1
      if (abs(a(j - 1, j)) < dprec) then
        p = j
        exit
      end if
    end do
  end subroutine find_nonzero_superdiagonal
 
  !!> @brief Checks if a matrix has a (near-)zero diagonal element.
  !!
  !<
  pure function has_zero_diagonal(A) result(has_zero)
    ! -- dummy
    real(dp), intent(in) :: a(:, :)
    logical(LGP) :: has_zero
    ! -- locals
    integer(I4B) :: m, n, i
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    has_zero = .false.
    do i = 1, min(m, n)
      if (abs(a(i, i)) < dprec) then
        has_zero = .true.
        exit
      end if
    end do
 
  end function has_zero_diagonal
 
  !!> @brief Reduces a rectangular matrix to square form using Givens rotations.
  !!
  !! This subroutine transforms a rectangular matrix `A` into a square matrix by applying Givens rotations
  !! to eliminate extra columns (when the number of columns exceeds the number of rows).
  !! The corresponding orthogonal transformations are accumulated in the matrix `Qt`.
  !<
  pure subroutine make_matrix_square(A, Qt)
    ! -- dummy
    real(dp), intent(inout), dimension(:, :) :: a
    real(dp), intent(inout), dimension(:, :), allocatable :: qt
    ! -- locals
    real(dp), dimension(2, 2) :: g
    integer(I4B) :: m, n, i
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    do i = m, 1, -1
      g = givens_rotation(a(i, i), a(i, m + 1))
      a(:, [i, m + 1]) = matmul(a(:, [i, m + 1]), transpose(g))
      qt([i, m + 1], :) = matmul(g, qt([i, m + 1], :))
    end do
 
  end subroutine make_matrix_square
 
  !!> @brief Sets small superdiagonal elements to zero for numerical stability.
  !!
  !! This subroutine scans the superdiagonal elements of the matrix `A` and sets to zero any element
  !! whose absolute value is less than or equal to a tolerance (relative to the neighboring diagonal elements).
  !! This helps to maintain numerical stability and ensures that the matrix is properly treated as bidiagonal
  !! in subsequent SVD steps.
  !!
  !! The criterion used is:
  !!   |A(j, j+1)| <= DPREC * (|A(j, j)| + |A(j+1, j+1)|)
  !! In such cases, A(j, j+1) is set to zero.
  !<
  pure subroutine clean_superdiagonal(A)
    ! -- dummy
    real(dp), intent(inout), dimension(:, :) :: a
    ! -- locals
    integer(I4B) :: n, m, j
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
 
    do j = 1, min(n, m) - 1
      if (abs(a(j, j + 1)) <= dprec * (abs(a(j, j)) + abs(a(j + 1, j + 1)))) then
        a(j, j + 1) = 0.0_dp
      end if
    end do
 
  end subroutine clean_superdiagonal
 
  !> @brief Singular Value Decomposition
  !!
  !! This method decomposes the matrix A into U, S and VT.
  !! It follows the algorithm as described by Golub and Reinsch.
  !!
  !! The first step is to decompose the matrix A into a bidiagonal matrix.
  !! This is done using Householder transformations.
  !! Then second step is to decompose the bidiagonal matrix into U, S and VT
  !! by repetitively applying the QR algorithm.
  !! If there is a zero on the diagonal or superdiagonal the matrix can be split
  !! into two smaller matrices and the QR algorithm can be applied to the smaller
  !! matrices.
  !!
  !! The matrix U is the eigenvectors of A*A^T
  !! The matrix VT is the eigenvectors of A^T*A
  !! The matrix S is the square root of the eigenvalues of A*A^T or A^T*A
  !!
  !<
  pure subroutine svd(A, U, S, VT)
    ! -- dummy
    real(dp), intent(in), dimension(:, :) :: a
    real(dp), intent(out), dimension(:, :), allocatable :: u
    real(dp), intent(out), dimension(:, :), allocatable :: s
    real(dp), intent(out), dimension(:, :), allocatable :: vt
    ! -- locals
    integer(I4B) :: i ! Loop counter for QR iterations
    integer(I4B) :: max_itr ! Maximum number of QR iterations allowed for convergence
    real(dp) :: error ! Convergence criterion (superdiagonal norm)
 
    integer(I4B) :: m ! Number of rows in input matrix A
    integer(I4B) :: n ! Number of columns in input matrix A
 
    integer(I4B) :: r ! Start index of active submatrix (for QR step)
    integer(I4B) :: q ! End index of active submatrix (for QR step)
 
    max_itr = 100
 
    m = size(a, dim=1) ! Number of rows
    n = size(a, dim=2) ! Number of columns
    s = a
 
    call bidiagonal_decomposition(s, u, vt)
 
    if (n > m) then
      call make_matrix_square(s, vt)
    end if
 
    do i = 1, max_itr
      call find_nonzero_superdiagonal(s, r, q)
 
      ! find zero diagonal
      if (has_zero_diagonal(s(r:q, r:q))) then
        ! write(*,*) 'Iteration: ', i, ' handle zero diagonal element'
        call handle_zero_diagonal(s(r:q, r:q), u(:, r:q), vt(r:q, :))
      else
        call bidiagonal_qr_decomposition(s(r:q, r:q), u(:, r:q), vt(r:q, :))
        call clean_superdiagonal(s(r:q, r:q))
 
        error = superdiagonal_norm(s)
        ! write(*,*) 'Iteration: ', i, ' Error: ', error
        if (error < dprec) exit
      end if
    end do
 
  end subroutine svd
 
end module