itkMathSVD.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkMathSVD_h
#define itkMathSVD_h
 
#include "itkMacro.h"
#include "itkMatrix.h"
#include "itkEigenDecompositionSignConvention.h"
#include "vnl/vnl_matrix.h"
#include "vnl/vnl_matrix_fixed.h"
#include "vnl/vnl_vector.h"
#include "vnl/vnl_vector_fixed.h"
 
#include "itk_eigen.h"
#include ITK_EIGEN(Dense)
 
#include <limits>
 
namespace itk
{
namespace Math
{
 
namespace detail
{
// rcond < 0 selects an automatic relative threshold of n * epsilon.
template <typename TReal>
TReal
ResolveRcond(TReal rcond, unsigned int n)
{
  return rcond < TReal{ 0 } ? static_cast<TReal>(n) * std::numeric_limits<TReal>::epsilon() : rcond;
}
 
// Moore-Penrose pseudo-inverse V diag(1/w) U^T, with singular values at or below
// rcond*max(w) treated as zero. Works for vnl_matrix and vnl_matrix_fixed.
template <typename TMatrix, typename TVector, typename TReal>
TMatrix
PseudoInverse(const TMatrix & U, const TVector & W, const TMatrix & V, TReal rcond)
{
  const unsigned int k = W.size();
  // W is sorted descending (Eigen guarantee), so W[0] is the max without an O(k) scan.
  const TReal tol = ResolveRcond(rcond, k) * W[0];
  TMatrix     scaledV = V;
  for (unsigned int col = 0; col < k; ++col)
  {
    const TReal s = (W[col] > tol) ? TReal{ 1 } / W[col] : TReal{ 0 };
    for (unsigned int i = 0; i < scaledV.rows(); ++i)
    {
      scaledV(i, col) *= s;
    }
  }
  return scaledV * U.transpose();
}
 
// Least-squares / minimum-norm solution x = V diag(1/w) U^T b of A x = b.
template <typename TMatrix, typename TVector, typename TReal>
TVector
SolveLinear(const TMatrix & U, const TVector & W, const TMatrix & V, const TVector & b, TReal rcond)
{
  const unsigned int n = W.size();
  const TReal        tol = ResolveRcond(rcond, n) * W[0]; // W sorted descending; W[0] == max
  TVector            utb = U.transpose() * b;
  for (unsigned int k = 0; k < n; ++k)
  {
    utb[k] = (W[k] > tol) ? utb[k] / W[k] : TReal{ 0 };
  }
  return V * utb;
}
 
template <typename TVector, typename TReal>
unsigned int
NumericalRank(const TVector & W, TReal rcond)
{
  const unsigned int n = W.size();
  const TReal        tol = ResolveRcond(rcond, n) * W[0]; // W sorted descending; W[0] == max
  unsigned int       count = 0;
  for (unsigned int k = 0; k < n; ++k)
  {
    if (W[k] > tol)
    {
      ++count;
    }
  }
  return count;
}
 
// Reconstruct U diag(w) V^T, treating singular values at or below rcond*max(w) as
// zero (a truncated, rank-reduced reconstruction of A).
template <typename TMatrix, typename TVector, typename TReal>
TMatrix
Recompose(const TMatrix & U, const TVector & W, const TMatrix & V, TReal rcond)
{
  const unsigned int k = W.size();
  const TReal        tol = ResolveRcond(rcond, k) * W[0]; // W sorted descending; W[0] == max
  TMatrix            scaledU = U;
  for (unsigned int col = 0; col < k; ++col)
  {
    const TReal s = (W[col] > tol) ? W[col] : TReal{ 0 };
    for (unsigned int i = 0; i < scaledU.rows(); ++i)
    {
      scaledU(i, col) *= s;
    }
  }
  return scaledU * V.transpose();
}
} // namespace detail

template <typename TReal, unsigned int VDim>
struct FixedSquareSVDResult
{
  vnl_matrix_fixed<TReal, VDim, VDim> U{};
  vnl_vector_fixed<TReal, VDim>       W{};
  vnl_matrix_fixed<TReal, VDim, VDim> V{};

  vnl_matrix_fixed<TReal, VDim, VDim>
  PseudoInverse(TReal rcond = TReal{ -1 }) const
  {
    return detail::PseudoInverse(U, W, V, rcond);
  }

  vnl_vector_fixed<TReal, VDim>
  Solve(const vnl_vector_fixed<TReal, VDim> & b, TReal rcond = TReal{ -1 }) const
  {
    return detail::SolveLinear(U, W, V, b, rcond);
  }

  unsigned int
  Rank(TReal rcond = TReal{ -1 }) const
  {
    return detail::NumericalRank(W, rcond);
  }

  vnl_matrix_fixed<TReal, VDim, VDim>
  Recompose(TReal rcond = TReal{ -1 }) const
  {
    return detail::Recompose(U, W, V, rcond);
  }
};

template <typename TReal>
struct SVDResult
{
  vnl_matrix<TReal> U{};
  vnl_vector<TReal> W{};
  vnl_matrix<TReal> V{};

  vnl_matrix<TReal>
  PseudoInverse(TReal rcond = TReal{ -1 }) const
  {
    return detail::PseudoInverse(U, W, V, rcond);
  }

  vnl_vector<TReal>
  Solve(const vnl_vector<TReal> & b, TReal rcond = TReal{ -1 }) const
  {
    return detail::SolveLinear(U, W, V, b, rcond);
  }

  unsigned int
  Rank(TReal rcond = TReal{ -1 }) const
  {
    return detail::NumericalRank(W, rcond);
  }

  vnl_matrix<TReal>
  Recompose(TReal rcond = TReal{ -1 }) const
  {
    return detail::Recompose(U, W, V, rcond);
  }
};
 
namespace detail
{
// Above this compile-time size the fixed overload delegates to the runtime path
// (keeps a large VDim off Eigen's stack-allocation limit).
constexpr unsigned int kFixedSVDMaxDim = 16;
 
// JacobiSVD+NoQRPreconditioner is fastest for small square inputs; BDCSVD is
// faster and more accurate for larger n.
constexpr unsigned int kJacobiMaxDim = 6;
 
// Runtime-sized square SVD: JacobiSVD + NoQRPreconditioner for small n, BDCSVD for
// larger n where Jacobi sweeps scale poorly. Throws on a failed decomposition.
template <typename TReal>
void
DynamicSquareSVDEigen(const TReal * inData, unsigned int n, TReal * uData, TReal * wData, TReal * vData)
{
  using RowMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
  using ColMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic>;
  constexpr int full = Eigen::ComputeFullU | Eigen::ComputeFullV;
 
  const Eigen::Map<const RowMajor> inMap(inData, n, n);
  Eigen::Map<RowMajor>             uMap(uData, n, n);
  Eigen::Map<RowMajor>             vMap(vData, n, n);
 
  const auto extract = [&](const auto & svd) {
    if (svd.info() != Eigen::Success)
    {
      itkGenericExceptionMacro("itk::Math::SVD failed; input is likely non-finite (NaN/Inf).");
    }
    uMap = svd.matrixU();
    vMap = svd.matrixV();
    for (unsigned int i = 0; i < n; ++i)
    {
      wData[i] = svd.singularValues()[i];
    }
  };
 
  if (n <= kJacobiMaxDim)
  {
    extract(Eigen::JacobiSVD < ColMajor, full | Eigen::NoQRPreconditioner > (inMap));
  }
  else
  {
    extract(Eigen::BDCSVD<ColMajor, full>(inMap));
  }
}
 
// NoQRPreconditioner skips the column-pivot QR (wasted for a square input); large
// VDim falls back to the runtime path. Throws on a failed/non-finite decomposition.
template <unsigned int VDim, typename TReal>
void
SquareSVDEigen(const TReal * inData, TReal * uData, TReal * wData, TReal * vData)
{
  if constexpr (VDim > kFixedSVDMaxDim)
  {
    DynamicSquareSVDEigen<TReal>(inData, VDim, uData, wData, vData);
  }
  else
  {
    // Compile-time sizes ≤ kFixedSVDMaxDim always take JacobiSVD: it stack-allocates,
    // fully unrolls, and beats BDCSVD's dynamic setup at these sizes (kJacobiMaxDim
    // gates only the runtime path, where BDCSVD's allocation pays off for larger n).
    using RowMajor = Eigen::Matrix<TReal, VDim, VDim, Eigen::RowMajor>;
    using ColMajor = Eigen::Matrix<TReal, VDim, VDim>;
    constexpr int options = Eigen::ComputeFullU | Eigen::ComputeFullV | Eigen::NoQRPreconditioner;
 
    const Eigen::Map<const RowMajor>          inMap(inData);
    const Eigen::JacobiSVD<ColMajor, options> svd(inMap);
    if (svd.info() != Eigen::Success)
    {
      itkGenericExceptionMacro("itk::Math::SVD failed; input is likely non-finite (NaN/Inf).");
    }
 
    Eigen::Map<RowMajor> uMap(uData);
    Eigen::Map<RowMajor> vMap(vData);
    uMap = svd.matrixU();
    vMap = svd.matrixV();
    for (unsigned int i = 0; i < VDim; ++i)
    {
      wData[i] = svd.singularValues()[i];
    }
  }
}
 
// Rectangular SVD via BDCSVD with thin U/V (Eigen's general engine for non-square
// inputs). For an m x n input: U is m x k, V is n x k, W has length k = min(m, n).
// Throws on a failed decomposition.
template <typename TReal>
void
RectangularSVDEigen(const TReal * inData,
                    unsigned int  rows,
                    unsigned int  cols,
                    TReal *       uData,
                    TReal *       wData,
                    TReal *       vData)
{
  using RowMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
  using ColMajor = Eigen::Matrix<TReal, Eigen::Dynamic, Eigen::Dynamic>;
  constexpr int      thin = Eigen::ComputeThinU | Eigen::ComputeThinV;
  const unsigned int k = (rows < cols) ? rows : cols;
 
  const Eigen::Map<const RowMajor>    inMap(inData, rows, cols);
  const Eigen::BDCSVD<ColMajor, thin> svd(inMap);
  if (svd.info() != Eigen::Success)
  {
    itkGenericExceptionMacro("itk::Math::SVD failed; input is likely non-finite (NaN/Inf).");
  }
  Eigen::Map<RowMajor> uMap(uData, rows, k);
  Eigen::Map<RowMajor> vMap(vData, cols, k);
  uMap = svd.matrixU();
  vMap = svd.matrixV();
  for (unsigned int i = 0; i < k; ++i)
  {
    wData[i] = svd.singularValues()[i];
  }
}
} // namespace detail

template <typename TReal, unsigned int VDim>
FixedSquareSVDResult<TReal, VDim>
SVD(const vnl_matrix_fixed<TReal, VDim, VDim> & A, bool canonicalizeSigns = true)
{
  FixedSquareSVDResult<TReal, VDim> result;
  detail::SquareSVDEigen<VDim>(A.data_block(), result.U.data_block(), result.W.data_block(), result.V.data_block());
  if (canonicalizeSigns)
  {
    itk::detail::CanonicalizeColumnSignsPaired(result.U, result.V);
  }
  return result;
}

template <typename TReal, unsigned int VDim>
FixedSquareSVDResult<TReal, VDim>
SVD(const Matrix<TReal, VDim, VDim> & A, bool canonicalizeSigns = true)
{
  return SVD<TReal, VDim>(A.GetVnlMatrix(), canonicalizeSigns);
}

template <typename TReal>
SVDResult<TReal>
SVD(const vnl_matrix<TReal> & A, bool canonicalizeSigns = true)
{
  const unsigned int rows = A.rows();
  const unsigned int cols = A.cols();
  if (rows == 0 || cols == 0)
  {
    itkGenericExceptionMacro("itk::Math::SVD requires a non-empty matrix.");
  }
  const unsigned int k = (rows < cols) ? rows : cols;
  SVDResult<TReal>   result;
  result.U.set_size(rows, k);
  result.V.set_size(cols, k);
  result.W.set_size(k);
  if (rows == cols)
  {
    detail::DynamicSquareSVDEigen<TReal>(
      A.data_block(), rows, result.U.data_block(), result.W.data_block(), result.V.data_block());
  }
  else
  {
    detail::RectangularSVDEigen<TReal>(
      A.data_block(), rows, cols, result.U.data_block(), result.W.data_block(), result.V.data_block());
  }
  if (canonicalizeSigns)
  {
    itk::detail::CanonicalizeColumnSignsPaired(result.U, result.V);
  }
  return result;
}
 
} // namespace Math
} // namespace itk
 
#endif // itkMathSVD_h