PseudoInverse.f90

Go to the documentation of this file.

module pseudoinversemodule
  use kindmodule, only: dp, i4b
  use constantsmodule, only: dprec
 
  use svdmodule, only: svd
 
  implicit none
  private
 
  public :: pinv
 
contains
 
  !!> @brief Computes the Moore-Penrose pseudoinverse of a matrix using SVD.
  !!
  !! This function computes the pseudoinverse \‍(A^+\‍) of the input matrix \‍(A\‍) using
  !! the singular value decomposition (SVD). Small singular values (less than `DPREC`)
  !! are set to zero to improve numerical stability. The pseudoinverse is useful for
  !! solving least-squares problems and for handling rank-deficient or non-square matrices.
  !!
  !! The result satisfies:
  !!     B = V * Sigma^+ * U^T
  !! where \‍(A = U \Sigma V^T\‍) is the SVD of \‍(A\‍), and \‍(\Sigma^+\‍) is the diagonal matrix
  !! with reciprocals of the nonzero singular values.
  !<
  function pinv(A) result(B)
    ! -- dummy
    real(dp), intent(in) :: a(:, :)
    real(dp) :: b(size(a, dim=2), size(a, dim=1)) !! The pseudoinverse of A (size n x m, where A is m x n)
    ! -- locals
    integer(I4B) :: m, n
    integer(I4B) :: pos
    real(dp), dimension(:, :), allocatable :: u
    real(dp), dimension(:, :), allocatable :: vt
    real(dp), dimension(:, :), allocatable :: sigma
 
    m = size(a, dim=1)
    n = size(a, dim=2)
 
    CALL svd(a, u, sigma, vt)
 
    ! Transform Sigma to its Sigma^+
    do pos = 1, min(m, n)
      if (dabs(sigma(pos, pos)) < dprec) then
        sigma(pos, pos) = 0.0_dp
      else
        sigma(pos, pos) = 1.0_dp / sigma(pos, pos)
      end if
    end do
 
    b = matmul(transpose(vt), matmul(transpose(sigma), transpose(u)))
 
  end function pinv
 
end module pseudoinversemodule