Class PsdMatrices

java.lang.Object
com.quantfinlib.risk.PsdMatrices

public final class PsdMatrices extends Object
Factorisation and repair of symmetric positive semi-definite matrices: the one place in the risk, simulation and copula code that decides what a "valid covariance" is.
  • choleskyPsd(double[][]) factors A = L·L' for a positive semi-definite A. A pivot that is zero to rounding (a rank-deficient factor covariance, two identical assets) gives a zero column, so the factor still reproduces A. A pivot that is clearly negative, or a zero pivot whose column is not zero, means the matrix is indefinite and the method throws. Tolerances are relative to the largest diagonal entry, so the result does not depend on units.
  • symmetricEigen(double[][]) is the cyclic Jacobi eigen-decomposition (Golub and Van Loan, 4th ed., Algorithm 8.5.3). It is O(n³) per sweep and suited to a few hundred dimensions at most.
  • nearestPsd(double[][]) clips negative eigenvalues to zero, which is the Frobenius-nearest PSD matrix (Higham 1988).
  • nearestCorrelation(double[][]) is Higham's (2002) alternating projections method with Dykstra's correction, unit weights: the Frobenius-nearest matrix that is PSD with unit diagonal.

References: N. J. Higham, "Computing the nearest correlation matrix - a problem from finance", IMA J. Numer. Anal. 22 (2002) 329-343; N. J. Higham, "Computing a nearest symmetric positive semidefinite matrix", Linear Algebra Appl. 103 (1988) 103-118.

  • Nested Class Summary

    Nested Classes
    Modifier and Type
    Class
    Description
    static final class 
    Eigenvalues (descending) and matching unit eigenvectors vectors[k].
  • Method Summary

    Modifier and Type
    Method
    Description
    static double[][]
    choleskyPsd(double[][] a)
    Lower-triangular L with A = L·L' for a symmetric PSD matrix.
    static double[][]
    nearestCorrelation(double[][] a)
    Nearest correlation matrix with tolerance 1e-12 and at most 10 000 iterations.
    static double[][]
    nearestCorrelation(double[][] a, double tol, int maxIter)
    Higham (2002) alternating projections with Dykstra's correction: R = Y − ΔS; X = P_S(R); ΔS = X − R; Y = P_U(X), where P_S clips eigenvalues at zero and P_U sets the diagonal to one.
    static double[][]
    nearestPsd(double[][] a)
    Frobenius-nearest PSD matrix: eigenvalues below zero set to zero.
    symmetricEigen(double[][] a)
    Cyclic Jacobi eigen-decomposition of a symmetric matrix.

    Methods inherited from class Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
  • Method Details

    • choleskyPsd

      public static double[][] choleskyPsd(double[][] a)
      Lower-triangular L with A = L·L' for a symmetric PSD matrix. See the class comment for how semi-definite and indefinite input are told apart.
      Throws:
      IllegalArgumentException - if a is not square, not symmetric, not finite, or indefinite
    • symmetricEigen

      public static PsdMatrices.Eigen symmetricEigen(double[][] a)
      Cyclic Jacobi eigen-decomposition of a symmetric matrix.
    • nearestPsd

      public static double[][] nearestPsd(double[][] a)
      Frobenius-nearest PSD matrix: eigenvalues below zero set to zero.
    • nearestCorrelation

      public static double[][] nearestCorrelation(double[][] a)
      Nearest correlation matrix with tolerance 1e-12 and at most 10 000 iterations.
    • nearestCorrelation

      public static double[][] nearestCorrelation(double[][] a, double tol, int maxIter)
      Higham (2002) alternating projections with Dykstra's correction: R = Y − ΔS; X = P_S(R); ΔS = X − R; Y = P_U(X), where P_S clips eigenvalues at zero and P_U sets the diagonal to one. Stops when ||Y − X||_F <= tol·||Y||_F.
      Throws:
      IllegalStateException - if it has not converged in maxIter