Matrix.java

  1. package org.drip.numerical.linearalgebra;

  3. /*
  4.  * -*- mode: java; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
  5.  */

  7. /*!
  8.  * Copyright (C) 2020 Lakshmi Krishnamurthy
  9.  * Copyright (C) 2019 Lakshmi Krishnamurthy
  10.  * Copyright (C) 2018 Lakshmi Krishnamurthy
  11.  * Copyright (C) 2017 Lakshmi Krishnamurthy
  12.  * Copyright (C) 2016 Lakshmi Krishnamurthy
  13.  * Copyright (C) 2015 Lakshmi Krishnamurthy
  14.  * Copyright (C) 2014 Lakshmi Krishnamurthy
  15.  * Copyright (C) 2013 Lakshmi Krishnamurthy
  16.  * 
  17.  *  This file is part of DROP, an open-source library targeting analytics/risk, transaction cost analytics,
  18.  *      asset liability management analytics, capital, exposure, and margin analytics, valuation adjustment
  19.  *      analytics, and portfolio construction analytics within and across fixed income, credit, commodity,
  20.  *      equity, FX, and structured products. It also includes auxiliary libraries for algorithm support,
  21.  *      numerical analysis, numerical optimization, spline builder, model validation, statistical learning,
  22.  *      and computational support.
  23.  *  
  24.  *      https://lakshmidrip.github.io/DROP/
  25.  *  
  26.  *  DROP is composed of three modules:
  27.  *  
  28.  *  - DROP Product Core - https://lakshmidrip.github.io/DROP-Product-Core/
  29.  *  - DROP Portfolio Core - https://lakshmidrip.github.io/DROP-Portfolio-Core/
  30.  *  - DROP Computational Core - https://lakshmidrip.github.io/DROP-Computational-Core/
  31.  * 
  32.  *  DROP Product Core implements libraries for the following:
  33.  *  - Fixed Income Analytics
  34.  *  - Loan Analytics
  35.  *  - Transaction Cost Analytics
  36.  * 
  37.  *  DROP Portfolio Core implements libraries for the following:
  38.  *  - Asset Allocation Analytics
  39.  *  - Asset Liability Management Analytics
  40.  *  - Capital Estimation Analytics
  41.  *  - Exposure Analytics
  42.  *  - Margin Analytics
  43.  *  - XVA Analytics
  44.  * 
  45.  *  DROP Computational Core implements libraries for the following:
  46.  *  - Algorithm Support
  47.  *  - Computation Support
  48.  *  - Function Analysis
  49.  *  - Model Validation
  50.  *  - Numerical Analysis
  51.  *  - Numerical Optimizer
  52.  *  - Spline Builder
  53.  *  - Statistical Learning
  54.  * 
  55.  *  Documentation for DROP is Spread Over:
  56.  * 
  57.  *  - Main                     => https://lakshmidrip.github.io/DROP/
  58.  *  - Wiki                     => https://github.com/lakshmiDRIP/DROP/wiki
  59.  *  - GitHub                   => https://github.com/lakshmiDRIP/DROP
  60.  *  - Repo Layout Taxonomy     => https://github.com/lakshmiDRIP/DROP/blob/master/Taxonomy.md
  61.  *  - Javadoc                  => https://lakshmidrip.github.io/DROP/Javadoc/index.html
  62.  *  - Technical Specifications => https://github.com/lakshmiDRIP/DROP/tree/master/Docs/Internal
  63.  *  - Release Versions         => https://lakshmidrip.github.io/DROP/version.html
  64.  *  - Community Credits        => https://lakshmidrip.github.io/DROP/credits.html
  65.  *  - Issues Catalog           => https://github.com/lakshmiDRIP/DROP/issues
  66.  *  - JUnit                    => https://lakshmidrip.github.io/DROP/junit/index.html
  67.  *  - Jacoco                   => https://lakshmidrip.github.io/DROP/jacoco/index.html
  68.  * 
  69.  *  Licensed under the Apache License, Version 2.0 (the "License");
  70.  *      you may not use this file except in compliance with the License.
  71.  *   
  72.  *  You may obtain a copy of the License at
  73.  *      http://www.apache.org/licenses/LICENSE-2.0
  74.  *  
  75.  *  Unless required by applicable law or agreed to in writing, software
  76.  *      distributed under the License is distributed on an "AS IS" BASIS,
  77.  *      WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  78.  *  
  79.  *  See the License for the specific language governing permissions and
  80.  *      limitations under the License.
  81.  */

  83. /**
  84.  * <i>Matrix</i> implements Matrix manipulation routines. It exports the following functionality:
  85.  *  <ul>
  86.  *      <li>
  87.  *          Matrix Inversion using Closed form solutions (for low-dimension matrices), or using Gaussian
  88.  *              elimination
  89.  *      </li>
  90.  *      <li>
  91.  *          Matrix Product
  92.  *      </li>
  93.  *      <li>
  94.  *          Matrix Diagonalization and Diagonal Pivoting
  95.  *      </li>
  96.  *      <li>
  97.  *          Matrix Regularization through Row Addition/Row Swap
  98.  *      </li>
  99.  *  </ul>
  100.  * 
  101.  * <br><br>
  102.  *  <ul>
  103.  *      <li><b>Module </b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/ComputationalCore.md">Computational Core Module</a></li>
  104.  *      <li><b>Library</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/NumericalAnalysisLibrary.md">Numerical Analysis Library</a></li>
  105.  *      <li><b>Project</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/numerical/README.md">Numerical Quadrature, Differentiation, Eigenization, Linear Algebra, and Utilities</a></li>
  106.  *      <li><b>Package</b> = <a href = "https://github.com/lakshmiDRIP/DROP/tree/master/src/main/java/org/drip/numerical/linearalgebra/README.md">Linear Algebra Matrix Transform Library</a></li>
  107.  *  </ul>
  108.  * <br><br>
  109.  *
  110.  * @author Lakshmi Krishnamurthy
  111.  */

  113. public class Matrix {

  115.     /**
  116.      * Lower Triangular Matrix
  117.      */

  119.     public static int LOWER_TRIANGULAR = 1;

  121.     /**
  122.      * Upper Triangular Matrix
  123.      */

  125.     public static int UPPER_TRIANGULAR = 2;

  127.     /**
  128.      * Lower+Upper Triangular Matrix
  129.      */

  131.     public static int LOWER_AND_UPPER_TRIANGULAR = 3;

  133.     /**
  134.      * Non Triangular Matrix
  135.      */

  137.     public static int NON_TRIANGULAR = 0;

  139.     /**
  140.      * Diagonalize the specified row in the source matrix, and apply comparable operations to the target
  141.      * 
  142.      * @param iQ Row in the Source Matrix
  143.      * @param aadblZ2XJack Source Matrix
  144.      * @param aadblZ2YJack Target Matrix
  145.      * 
  146.      * @return TRUE - Diagonalization was successful
  147.      */

  149.     public static final boolean DiagonalizeRow (
  150.         final int iQ,
  151.         final double[][] aadblZ2XJack,
  152.         final double[][] aadblZ2YJack)
  153.     {
  154.         if (0. != aadblZ2XJack[iQ][iQ]) return true;

  156.         int iSize = aadblZ2XJack.length;
  157.         int iP = iSize - 1;

  159.         while (0. == aadblZ2XJack[iP][iQ] && iP >= 0) --iP;

  161.         if (0 > iP) return false;

  163.         for (int j = 0; j < iSize; ++j)
  164.             aadblZ2XJack[iQ][j] += aadblZ2XJack[iP][j];

  166.         aadblZ2YJack[iQ][iP] += 1.;
  167.         return true;
  168.     }

  170.     /**
  171.      * Compute the Product of an Input Matrix and a Column
  172.      * 
  173.      * @param aadblA Matrix A
  174.      * @param adblB Array B
  175.      * 
  176.      * @return The Product
  177.      */

  179.     public static final double[] Product (
  180.         final double[][] aadblA,
  181.         final double[] adblB)
  182.     {
  183.         if (null == aadblA || null == adblB) return null;

  185.         int iNumACol = aadblA[0].length;
  186.         int iNumProductCol = adblB.length;
  187.         int iNumProductRow = aadblA.length;
  188.         double[] adblProduct = new double[iNumProductRow];

  190.         if (0 == iNumACol || iNumACol != adblB.length || 0 == iNumProductRow || 0 == iNumProductCol)
  191.             return null;

  193.         for (int iRow = 0; iRow < iNumProductRow; ++iRow) {
  194.             adblProduct[iRow] = 0.;

  196.             for (int i = 0; i < iNumACol; ++i) {
  197.                 if (!org.drip.numerical.common.NumberUtil.IsValid (aadblA[iRow][i]) ||
  198.                     !org.drip.numerical.common.NumberUtil.IsValid (adblB[i]))
  199.                     return null;

  201.                 adblProduct[iRow] += aadblA[iRow][i] * adblB[i];
  202.             }
  203.         }

  205.         return adblProduct;
  206.     }

  208.     /**
  209.      * Compute the Product of an input column and a matrix
  210.      * 
  211.      * @param adblA Column A
  212.      * @param aadblB Matrix B
  213.      * 
  214.      * @return The Product
  215.      */

  217.     public static final double[][] Product (
  218.         final double[] adblA,
  219.         final double[][] aadblB)
  220.     {
  221.         if (null == adblA || null == aadblB) return null;

  223.         int iNumACol = adblA.length;
  224.         int iNumProductCol = aadblB.length;
  225.         double[][] aadblProduct = new double[iNumACol][iNumProductCol];

  227.         if (0 == iNumACol || iNumACol != aadblB.length || 0 == iNumProductCol) return null;

  229.         for (int iRow = 0; iRow < iNumACol; ++iRow) {
  230.             for (int iCol = 0; iCol < iNumProductCol; ++iCol) {
  231.                 aadblProduct[iRow][iCol] = 0.;

  233.                 for (int i = 0; i < iNumACol; ++i) {
  234.                     if (!org.drip.numerical.common.NumberUtil.IsValid (adblA[iRow]) ||
  235.                         !org.drip.numerical.common.NumberUtil.IsValid (aadblB[i][iCol]))
  236.                         return null;

  238.                     aadblProduct[iRow][iCol] += adblA[iRow] * aadblB[i][iCol];
  239.                 }
  240.             }
  241.         }

  243.         return aadblProduct;
  244.     }

  246.     /**
  247.      * Compute the Product of the input matrices
  248.      * 
  249.      * @param aadblA Matrix A
  250.      * @param aadblB Matrix B
  251.      * 
  252.      * @return The Product
  253.      */

  255.     public static final double[][] Product (
  256.         final double[][] aadblA,
  257.         final double[][] aadblB)
  258.     {
  259.         if (null == aadblA || null == aadblB) return null;

  261.         int iNumACol = aadblA[0].length;
  262.         int iNumProductRow = aadblA.length;
  263.         int iNumProductCol = aadblB[0].length;
  264.         double[][] aadblProduct = new double[iNumProductRow][iNumProductCol];

  266.         if (0 == iNumACol || iNumACol != aadblB.length || 0 == iNumProductRow || 0 == iNumProductCol)
  267.             return null;

  269.         for (int iRow = 0; iRow < iNumProductRow; ++iRow) {
  270.             for (int iCol = 0; iCol < iNumProductCol; ++iCol) {
  271.                 aadblProduct[iRow][iCol] = 0.;

  273.                 for (int i = 0; i < iNumACol; ++i) {
  274.                     if (!org.drip.numerical.common.NumberUtil.IsValid (aadblA[iRow][i]) ||
  275.                         !org.drip.numerical.common.NumberUtil.IsValid (aadblB[i][iCol]))
  276.                         return null;

  278.                     aadblProduct[iRow][iCol] += aadblA[iRow][i] * aadblB[i][iCol];
  279.                 }
  280.             }
  281.         }

  283.         return aadblProduct;
  284.     }

  286.     /**
  287.      * Make a Square Diagonal Matrix from a Row
  288.      * 
  289.      * @param adblA The Row Array
  290.      * 
  291.      * @return The corresponding Square Diagonal Matrix
  292.      */

  294.     public static final double[][] MakeSquareDiagonal (
  295.         final double[] adblA)
  296.     {
  297.         if (null == adblA) return null;

  299.         int iNumElement = adblA.length;
  300.         double[][] aadblDiagonal = 0 == iNumElement ? null : new double[iNumElement][iNumElement];

  302.         if (0 == iNumElement) return null;

  304.         for (int i = 0; i < iNumElement; ++i) {
  305.             for (int j = 0; j < iNumElement; ++j)
  306.                 aadblDiagonal[i][j] = i == j ? adblA[i] : 0.;
  307.         }

  309.         return aadblDiagonal;
  310.     }

  312.     /**
  313.      * Invert a 2D Matrix using Cramer's Rule
  314.      * 
  315.      * @param aadblA Input 2D Matrix
  316.      * 
  317.      * @return The Inverted Matrix
  318.      */

  320.     public static final double[][] Invert2DMatrixUsingCramerRule (
  321.         final double[][] aadblA)
  322.     {
  323.         if (null == aadblA || 2 != aadblA.length || 2 != aadblA[0].length) return null;

  325.         for (int i = 0; i < 2; ++i) {
  326.             for (int j = 0; j < 2; ++j) {
  327.                 if (!org.drip.numerical.common.NumberUtil.IsValid (aadblA[i][j])) return null;
  328.             }
  329.         }

  331.         double dblScale = aadblA[0][0] * aadblA[1][1] - aadblA[0][1] * aadblA[1][0];

  333.         if (0. == dblScale) return null;

  335.         return new double[][] {{aadblA[1][1] / dblScale, -1. * aadblA[0][1] / dblScale}, {-1. * aadblA[1][0]
  336.             / dblScale, aadblA[0][0] / dblScale}};
  337.     }

  339.     /**
  340.      * Regularize the specified diagonal entry of the input matrix using Row Swapping
  341.      * 
  342.      * @param mct The Input Matrix Complement Transform
  343.      * 
  344.      * @return The Regularization was successful
  345.      */

  347.     public static final boolean RegularizeUsingRowSwap (
  348.         final org.drip.numerical.linearalgebra.MatrixComplementTransform mct)
  349.     {
  350.         if (null == mct) return false;

  352.         int iSize = mct.size();

  354.         double[][] aadblSource = mct.getSource();

  356.         double[][] aadblComplement = mct.getComplement();

  358.         for (int iDiagonal = 0; iDiagonal < iSize; ++iDiagonal) {
  359.             if (0. == aadblSource[iDiagonal][iDiagonal]) {
  360.                 int iSwapRow = iSize - 1;

  362.                 while (iSwapRow >= 0 && (0. == aadblSource[iSwapRow][iDiagonal] || 0. ==
  363.                     aadblSource[iDiagonal][iSwapRow]))
  364.                     --iSwapRow;

  366.                 if (0 > iSwapRow) {
  367.                     iSwapRow = 0;

  369.                     while (iSwapRow < iSize && 0. == aadblSource[iSwapRow][iDiagonal])
  370.                         ++iSwapRow;

  372.                     if (iSwapRow >= iSize) return false;
  373.                 }

  375.                 for (int iCol = 0; iCol < iSize; ++iCol) {
  376.                     double dblComplementDiagonalEntry = aadblComplement[iDiagonal][iCol];
  377.                     aadblComplement[iDiagonal][iCol] = aadblComplement[iSwapRow][iCol];
  378.                     aadblComplement[iSwapRow][iCol] = dblComplementDiagonalEntry;
  379.                     double dblSourceDiagonalEntry = aadblSource[iDiagonal][iCol];
  380.                     aadblSource[iDiagonal][iCol] = aadblSource[iSwapRow][iCol];
  381.                     aadblSource[iSwapRow][iCol] = dblSourceDiagonalEntry;
  382.                 }
  383.             }
  384.         }

  386.         /* for (int iDiagonal = 0; iDiagonal < iSize; ++iDiagonal) {
  387.             if (0. == aadblSource[iDiagonal][iDiagonal]) {
  388.                 org.drip.quant.common.NumberUtil.Print2DArray ("ZERO DIAG!", aadblSource, false);

  390.                 return false;
  391.             }
  392.         } */

  394.         return true;
  395.     }

  397.     /**
  398.      * Regularize the specified diagonal entry of the input matrix using Row Addition
  399.      * 
  400.      * @param mct The Input Matrix Complement Transform
  401.      * 
  402.      * @return The Regularization was successful
  403.      */

  405.     public static final boolean RegularizeUsingRowAddition (
  406.         final org.drip.numerical.linearalgebra.MatrixComplementTransform mct)
  407.     {
  408.         if (null == mct) return false;

  410.         int iSize = mct.size();

  412.         double[][] aadblSource = mct.getSource();

  414.         double[][] aadblComplement = mct.getComplement();

  416.         for (int iDiagonal = 0; iDiagonal < iSize; ++iDiagonal) {
  417.             if (0. == aadblSource[iDiagonal][iDiagonal]) {
  418.                 int iPivotRow = iSize - 1;

  420.                 while (0. == aadblSource[iPivotRow][iDiagonal] && iPivotRow >= 0) --iPivotRow;

  422.                 if (0 > iPivotRow) return false;

  424.                 for (int iCol = 0; iCol < iSize; ++iCol) {
  425.                     aadblSource[iDiagonal][iCol] += aadblSource[iPivotRow][iCol];
  426.                     aadblComplement[iDiagonal][iCol] += aadblComplement[iPivotRow][iCol];
  427.                 }
  428.             }
  429.         }

  431.         return true;
  432.     }

  434.     /**
  435.      * Pivot the Diagonal of the Input Matrix
  436.      * 
  437.      * @param aadblA The Input Matrix
  438.      * 
  439.      * @return The Matrix Complement Transform Instance
  440.      */

  442.     public static final org.drip.numerical.linearalgebra.MatrixComplementTransform PivotDiagonal (
  443.         final double[][] aadblA)
  444.     {
  445.         if (null == aadblA) return null;

  447.         int iSize = aadblA.length;
  448.         double[][] aadblSource = new double[iSize][iSize];
  449.         double[][] aadblComplement = new double[iSize][iSize];
  450.         org.drip.numerical.linearalgebra.MatrixComplementTransform mctOut = null;

  452.         if (0 == iSize || null == aadblA[0] || iSize != aadblA[0].length) return null;

  454.         for (int i = 0; i < iSize; ++i) {
  455.             for (int j = 0; j < iSize; ++j) {
  456.                 aadblSource[i][j] = aadblA[i][j];
  457.                 aadblComplement[i][j] = i == j ? 1. : 0.;
  458.             }
  459.         }

  461.         try {
  462.             mctOut = new org.drip.numerical.linearalgebra.MatrixComplementTransform (aadblSource,
  463.                 aadblComplement);
  464.         } catch (java.lang.Exception e) {
  465.             e.printStackTrace();

  467.             return null;
  468.         }

  470.         return RegularizeUsingRowSwap (mctOut) ? mctOut : null;
  471.     }

  473.     /**
  474.      * Invert the Source Matrix using Gaussian Elimination
  475.      * 
  476.      * @param aadblSource Source Matrix
  477.      * 
  478.      * @return The Inverted Matrix
  479.      */

  481.     public static final double[][] InvertUsingGaussianElimination (
  482.         final double[][] aadblSource)
  483.     {
  484.         org.drip.numerical.linearalgebra.MatrixComplementTransform mctRegularized =
  485.             org.drip.numerical.linearalgebra.Matrix.PivotDiagonal (aadblSource);

  487.         if (null == mctRegularized) return null;

  489.         double[][] aadblRegularizedSource = mctRegularized.getSource();

  491.         double[][] aadblRegularizedInverse = mctRegularized.getComplement();

  493.         int iSize = aadblRegularizedSource.length;

  495.         for (int iDiagonal = 0; iDiagonal < iSize; ++iDiagonal) {
  496.             if (0. == aadblRegularizedSource[iDiagonal][iDiagonal]) return null;

  498.             for (int iRow = 0; iRow < iSize; ++iRow) {
  499.                 if (iRow == iDiagonal || 0. == aadblRegularizedSource[iRow][iDiagonal]) continue;

  501.                 double dblColEntryEliminatorRatio = aadblRegularizedSource[iDiagonal][iDiagonal] /
  502.                     aadblRegularizedSource[iRow][iDiagonal];

  504.                 for (int iCol = 0; iCol < iSize; ++iCol) {
  505.                     aadblRegularizedSource[iRow][iCol] = aadblRegularizedSource[iRow][iCol] *
  506.                         dblColEntryEliminatorRatio - aadblRegularizedSource[iDiagonal][iCol];
  507.                     aadblRegularizedInverse[iRow][iCol] = aadblRegularizedInverse[iRow][iCol] *
  508.                         dblColEntryEliminatorRatio - aadblRegularizedInverse[iDiagonal][iCol];
  509.                 }
  510.             }
  511.         }

  513.         for (int iDiagonal = 0; iDiagonal < iSize; ++iDiagonal) {
  514.             double dblDiagScaleDown = aadblRegularizedSource[iDiagonal][iDiagonal];

  516.             if (0. == dblDiagScaleDown) return null;

  518.             for (int iCol = 0; iCol < iSize; ++iCol) {
  519.                 aadblRegularizedSource[iDiagonal][iCol] /= dblDiagScaleDown;
  520.                 aadblRegularizedInverse[iDiagonal][iCol] /= dblDiagScaleDown;
  521.             }
  522.         }

  524.         return aadblRegularizedInverse;
  525.     }

  527.     /**
  528.      * Invert the input matrix using the specified Method
  529.      * 
  530.      * @param aadblA Input Matrix
  531.      * @param strMethod The Inversion Method
  532.      * 
  533.      * @return The Inverted Matrix
  534.      */

  536.     public static final double[][] Invert (
  537.         final double[][] aadblA,
  538.         final java.lang.String strMethod)
  539.     {
  540.         if (null == aadblA) return null;

  542.         int iSize = aadblA.length;
  543.         double[][] aadblAInv = null;
  544.         double[][] aadblASource = new double[iSize][iSize];
  545.         double[][] aadblZ2YJack = new double[iSize][iSize];

  547.         if (0 == iSize || iSize != aadblA[0].length) return null;

  549.         for (int i = 0; i < iSize; ++i) {
  550.             for (int j = 0; j < iSize; ++j) {
  551.                 if (!org.drip.numerical.common.NumberUtil.IsValid (aadblASource[i][j] = aadblA[i][j]))
  552.                     return null;

  554.                 aadblZ2YJack[i][j] = i == j ? 1. : 0.;
  555.             }
  556.         }

  558.         for (int i = 0; i < iSize; ++i) {
  559.             if (0. == aadblASource[i][i] && !DiagonalizeRow (i, aadblASource, aadblZ2YJack)) return null;
  560.         }

  562.         if (null == strMethod || strMethod.isEmpty() || strMethod.equalsIgnoreCase ("GaussianElimination"))
  563.             aadblAInv = InvertUsingGaussianElimination (aadblASource);

  565.         if (null == aadblAInv || iSize != aadblAInv.length || iSize != aadblAInv[0].length) return null;

  567.         return Product (aadblAInv, aadblZ2YJack);
  568.     }

  570.     /**
  571.      * Compute the Rank of the Matrix
  572.      * 
  573.      * @param aadblSource Source Matrix
  574.      * 
  575.      * @return The Rank of the Matrix
  576.      * 
  577.      * @throws java.lang.Exception Thrown if the Rank Cannot be computed
  578.      */

  580.     public static final int Rank (
  581.         final double[][] aadblSource)
  582.         throws java.lang.Exception
  583.     {
  584.         if (null == aadblSource) return 0;

  586.         int iNumRow = aadblSource.length;

  588.         if (iNumRow == 0) return 0;

  590.         int iNumCol = aadblSource[0].length;

  592.         for (int iScanRow = 0; iScanRow < iNumRow; ++iScanRow) {
  593.             if (!org.drip.numerical.common.NumberUtil.IsValid (aadblSource[iScanRow]))
  594.                 throw new java.lang.Exception ("Matrix::Rank => Invalid Inputs");
  595.         }

  597.         double[][] aadblRegularizedSource = iNumRow < iNumCol ?
  598.             org.drip.numerical.linearalgebra.Matrix.Transpose (aadblSource) : aadblSource;

  600.         int iNumDependentRow = 0;
  601.         int iProcessedRow = aadblRegularizedSource.length;
  602.         int iProcessedCol = aadblRegularizedSource[0].length;

  604.         if (1 == iNumRow || 1 == iNumCol) return iProcessedRow;

  606.         for (int iScanRow = 0; iScanRow < iProcessedCol; ++iScanRow) {
  607.             for (int iRow = 0; iRow < iProcessedCol; ++iRow) {
  608.                 if (iRow == iScanRow || 0. == aadblRegularizedSource[iRow][iScanRow]) continue;

  610.                 double dblColEntryEliminatorRatio = aadblRegularizedSource[iScanRow][iScanRow] /
  611.                     aadblRegularizedSource[iRow][iScanRow];

  613.                 for (int iCol = 0; iCol < iProcessedCol; ++iCol)
  614.                     aadblRegularizedSource[iRow][iCol] = aadblRegularizedSource[iRow][iCol] *
  615.                         dblColEntryEliminatorRatio - aadblRegularizedSource[iScanRow][iCol];
  616.             }
  617.         }

  619.         for (int iScanRow = 0; iScanRow < iProcessedCol; ++iScanRow) {
  620.             if (0. == org.drip.numerical.linearalgebra.Matrix.Modulus (aadblRegularizedSource[iScanRow]))
  621.                 ++iNumDependentRow;
  622.         }

  624.         return iProcessedRow - iNumDependentRow;
  625.     }

  627.     /**
  628.      * Transpose the specified Square Matrix
  629.      * 
  630.      * @param aadblA The Input Square Matrix
  631.      * 
  632.      * @return The Transpose of the Square Matrix
  633.      */

  635.     public static final double[][] Transpose (
  636.         final double[][] aadblA)
  637.     {
  638.         if (null == aadblA) return null;

  640.         int iRowSize = aadblA.length;

  642.         if (0 == iRowSize || null == aadblA[0]) return null;

  644.         int iColSize = aadblA[0].length;
  645.         double[][] aadblATranspose = new double[iColSize][iRowSize];

  647.         if (0 == iColSize) return null;

  649.         for (int i = 0; i < iColSize; ++i) {
  650.             for (int j = 0; j < iRowSize; ++j)
  651.                 aadblATranspose[i][j] = aadblA[j][i];
  652.         }

  654.         return aadblATranspose;
  655.     }

  657.     /**
  658.      * Compute the Cholesky-Banachiewicz Factorization of the specified Matrix.
  659.      * 
  660.      * @param aadblA The Input Matrix
  661.      * 
  662.      * @return The Factorized Matrix
  663.      */

  665.     public static final double[][] CholeskyBanachiewiczFactorization (
  666.         final double[][] aadblA)
  667.     {
  668.         if (null == aadblA) return null;

  670.         int iSize = aadblA.length;
  671.         double[][] aadblL = new double[iSize][iSize];

  673.         if (0 == iSize || null == aadblA[0] || iSize != aadblA[0].length) return null;

  675.         for (int i = 0; i < iSize; ++i) {
  676.             for (int j = 0; j < iSize; ++j) {
  677.                 aadblL[i][j] = 0.;

  679.                 if (i == j) {
  680.                     for (int k = 0; k < j; ++k)
  681.                         aadblL[j][j] -= aadblL[j][k] * aadblL[j][k];

  683.                     aadblL[j][j] = java.lang.Math.sqrt (aadblL[j][j] + aadblA[j][j]);
  684.                 } else if (i > j) {
  685.                     for (int k = 0; k < j; ++k)
  686.                         aadblL[i][j] -= aadblL[i][k] * aadblL[j][k];

  688.                     aadblL[i][j] = (aadblA[i][j] + aadblL[i][j]) / aadblL[j][j];
  689.                 }
  690.             }
  691.         }

  693.         return aadblL;
  694.     }

  696.     /**
  697.      * Dot Product of Vectors A and E
  698.      * 
  699.      * @param adblA Vector A
  700.      * @param adblE Vector E
  701.      * 
  702.      * @return The Dot Product
  703.      * 
  704.      * @throws java.lang.Exception Thrown if the Dot-Product cannot be computed
  705.      */

  707.     public static final double DotProduct (
  708.         final double[] adblA,
  709.         final double[] adblE)
  710.         throws java.lang.Exception
  711.     {
  712.         if (null == adblA || null == adblE)
  713.             throw new java.lang.Exception ("Matrix::DotProduct => Invalid Inputs!");

  715.         int iSize = adblA.length;
  716.         double dblDotProduct = 0.;

  718.         if (0 == iSize || iSize != adblE.length)
  719.             throw new java.lang.Exception ("Matrix::DotProduct => Invalid Inputs!");

  721.         for (int i = 0; i < iSize; ++i)
  722.             dblDotProduct += adblE[i] * adblA[i];

  724.         return dblDotProduct;
  725.     }

  727.     /**
  728.      * Compute the Cross Product between the Specified Vectors
  729.      * 
  730.      * @param vector1 Vector #1
  731.      * @param vector2 Vector #2
  732.      * 
  733.      * @return The Cross Product
  734.      */

  736.     public static final double[][] CrossProduct (
  737.         final double[] vector1,
  738.         final double[] vector2)
  739.     {
  740.         if (null == vector1 || null == vector2)
  741.         {
  742.             return null;
  743.         }

  745.         int size1 = vector1.length;
  746.         int size2 = vector2.length;
  747.         double[][] crossProduct = 0 == size1 || 0 == size2 ? null : new double[size1][size2];

  749.         if (null == crossProduct)
  750.         {
  751.             return null;
  752.         }

  754.         for (int index1 = 0; index1 < size1; ++index1)
  755.         {
  756.             for (int index2 = 0; index2 < size2; ++index2)
  757.             {
  758.                 crossProduct[index1][index2] = vector1[index1] * vector2[index2];
  759.             }
  760.         }

  762.         return crossProduct;
  763.     }

  765.     /**
  766.      * Project the Vector A along the Vector E
  767.      * 
  768.      * @param adblA Vector A
  769.      * @param adblE Vector E
  770.      * 
  771.      * @return The Vector of Projection of A along E
  772.      */

  774.     public static final double[] Project (
  775.         final double[] adblA,
  776.         final double[] adblE)
  777.     {
  778.         if (null == adblA || null == adblE) return null;

  780.         int iSize = adblA.length;
  781.         double dblProjection = java.lang.Double.NaN;
  782.         double[] adblProjectAOnE = new double[iSize];

  784.         if (0 == iSize || iSize != adblE.length) return null;

  786.         try {
  787.             dblProjection = DotProduct (adblA, adblE) / DotProduct (adblE, adblE);
  788.         } catch (java.lang.Exception e) {
  789.             e.printStackTrace();

  791.             return null;
  792.         }

  794.         for (int i = 0; i < iSize; ++i)
  795.             adblProjectAOnE[i] = adblE[i] * dblProjection;

  797.         return adblProjectAOnE;
  798.     }

  800.     /**
  801.      * Compute the Sum of the Input Vector
  802.      * 
  803.      * @param adbl The Input Vector
  804.      * 
  805.      * @return TRUE - The Sum of the Input Vector
  806.      * 
  807.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  808.      */

  810.     public static final double Sum (
  811.         final double[] adbl)
  812.         throws java.lang.Exception
  813.     {
  814.         if (null == adbl || !org.drip.numerical.common.NumberUtil.IsValid (adbl))
  815.             throw new java.lang.Exception ("Matrix::Sum => Invalid Inputs");

  817.         double dblSum = 0.;
  818.         int iSize = adbl.length;

  820.         if (0 == iSize) throw new java.lang.Exception ("Matrix::Sum => Invalid Inputs");

  822.         for (int i = 0; i < iSize; ++i)
  823.             dblSum += adbl[i];

  825.         return dblSum;
  826.     }

  828.     /**
  829.      * Compute the Modulus of the Input Vector
  830.      * 
  831.      * @param adbl The Input Vector
  832.      * 
  833.      * @return TRUE - The Modulus of the Input Vector
  834.      * 
  835.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  836.      */

  838.     public static final double Modulus (
  839.         final double[] adbl)
  840.         throws java.lang.Exception
  841.     {
  842.         if (null == adbl || !org.drip.numerical.common.NumberUtil.IsValid (adbl))
  843.             throw new java.lang.Exception ("Matrix::Modulus => Invalid Inputs");

  845.         double dblModulus = 0.;
  846.         int iSize = adbl.length;

  848.         if (0 == iSize) throw new java.lang.Exception ("Matrix::Modulus => Invalid Inputs");

  850.         for (int i = 0; i < iSize; ++i)
  851.             dblModulus += adbl[i] * adbl[i];

  853.         return java.lang.Math.sqrt (dblModulus);
  854.     }

  856.     /**
  857.      * Indicate if the Array Entries are Positive or Zero
  858.      * 
  859.      * @param adbl The Array
  860.      * 
  861.      * @return TRUE - The Array Entries are Positive or Zero
  862.      * 
  863.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  864.      */

  866.     public static final boolean PositiveOrZero (
  867.         final double[] adbl)
  868.         throws java.lang.Exception
  869.     {
  870.         if (null == adbl || !org.drip.numerical.common.NumberUtil.IsValid (adbl))
  871.             throw new java.lang.Exception ("Matrix::PositiveOrZero => Invalid Inputs");

  873.         int iSize = adbl.length;

  875.         if (0 == iSize) throw new java.lang.Exception ("Matrix::PositiveOrZero => Invalid Inputs");

  877.         for (int i = 0; i < iSize; ++i) {
  878.             if (0. > adbl[i]) return false;
  879.         }

  881.         return true;
  882.     }

  884.     /**
  885.      * Indicate if the Array Entries are Negative or Zero
  886.      * 
  887.      * @param adbl The Array
  888.      * 
  889.      * @return The Array Entries are Negative or Zero
  890.      * 
  891.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  892.      */

  894.     public static final boolean NegativeOrZero (
  895.         final double[] adbl)
  896.         throws java.lang.Exception
  897.     {
  898.         if (null == adbl || !org.drip.numerical.common.NumberUtil.IsValid (adbl))
  899.             throw new java.lang.Exception ("Matrix::NegativeOrZero => Invalid Inputs");

  901.         int iSize = adbl.length;

  903.         if (0 == iSize)  throw new java.lang.Exception ("Matrix::NegativeOrZero => Invalid Inputs");

  905.         for (int i = 0; i < iSize; ++i) {
  906.             if (0. < adbl[i]) return false;
  907.         }

  909.         return true;
  910.     }

  912.     /**
  913.      * Indicate if the Array Entries are Positive Linearly Independent
  914.      * 
  915.      * @param adbl The Array
  916.      * 
  917.      * @return TRUE - The Array Entries are Positive Linearly Independent
  918.      * 
  919.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  920.      */

  922.     public static final boolean PositiveLinearlyIndependent (
  923.         final double[] adbl)
  924.         throws java.lang.Exception
  925.     {
  926.         return !PositiveOrZero (adbl) && !NegativeOrZero (adbl);
  927.     }

  929.     /**
  930.      * Normalize the Input Vector
  931.      * 
  932.      * @param adbl The Input Vector
  933.      * 
  934.      * @return The Normalized Vector
  935.      */

  937.     public static final double[] Normalize (
  938.         final double[] adbl)
  939.     {
  940.         if (null == adbl) return null;

  942.         double dblNorm = 0.;
  943.         int iSize = adbl.length;
  944.         double[] adblNormalized = new double[iSize];

  946.         if (0 == iSize) return null;

  948.         for (int i = 0; i < iSize; ++i)
  949.             dblNorm += adbl[i] * adbl[i];

  951.         dblNorm = java.lang.Math.sqrt (dblNorm);

  953.         for (int i = 0; i < iSize; ++i)
  954.             adblNormalized[i] = adbl[i] / dblNorm;

  956.         return adblNormalized;
  957.     }

  959.     /**
  960.      * Orthogonalize the Specified Matrix Using the Graham-Schmidt Method
  961.      * 
  962.      * @param aadblV The Input Matrix
  963.      * 
  964.      * @return The Orthogonalized Matrix
  965.      */

  967.     public static final double[][] GrahamSchmidtOrthogonalization (
  968.         final double[][] aadblV)
  969.     {
  970.         if (null == aadblV) return null;

  972.         int iSize = aadblV.length;
  973.         double[][] aadblU = new double[iSize][iSize];

  975.         if (0 == iSize || null == aadblV[0] || iSize != aadblV[0].length) return null;

  977.         for (int i = 0; i < iSize; ++i) {
  978.             for (int j = 0; j < iSize; ++j)
  979.                 aadblU[i][j] = aadblV[i][j];

  981.             for (int j = 0; j < i; ++j) {
  982.                 double dblProjectionAmplitude = java.lang.Double.NaN;

  984.                 try {
  985.                     dblProjectionAmplitude = DotProduct (aadblV[i], aadblU[j]) / DotProduct (aadblU[j],
  986.                         aadblU[j]);
  987.                 } catch (java.lang.Exception e) {
  988.                     e.printStackTrace();

  990.                     return null;
  991.                 }

  993.                 for (int k = 0; k < iSize; ++k)
  994.                     aadblU[i][k] -= dblProjectionAmplitude * aadblU[j][k];
  995.             }
  996.         }

  998.         return aadblU;
  999.     }

  1001.     /**
  1002.      * Orthonormalize the Specified Matrix Using the Graham-Schmidt Method
  1003.      * 
  1004.      * @param aadblV The Input Matrix
  1005.      * 
  1006.      * @return The Orthonormalized Matrix
  1007.      */

  1009.     public static final double[][] GrahamSchmidtOrthonormalization (
  1010.         final double[][] aadblV)
  1011.     {
  1012.         double[][] aadblVOrthogonal = GrahamSchmidtOrthogonalization (aadblV);

  1014.         if (null == aadblVOrthogonal) return null;

  1016.         int iSize = aadblVOrthogonal.length;

  1018.         double[][] aadblVOrthonormal = new double[iSize][];

  1020.         for (int i = 0; i < iSize; ++i)
  1021.             aadblVOrthonormal[i] = Normalize (aadblVOrthogonal[i]);

  1023.         return aadblVOrthonormal;
  1024.     }

  1026.     /**
  1027.      * Perform a QR Decomposition on the Input Matrix
  1028.      * 
  1029.      * @param aadblA The Input Matrix
  1030.      * 
  1031.      * @return The Output of QR Decomposition
  1032.      */

  1034.     public static final org.drip.numerical.linearalgebra.QR QRDecomposition (
  1035.         final double[][] aadblA)
  1036.     {
  1037.         double[][] aadblQ = GrahamSchmidtOrthonormalization (aadblA);

  1039.         if (null == aadblQ) return null;

  1041.         int iSize = aadblQ.length;
  1042.         double[][] aadblR = new double[iSize][iSize];

  1044.         try {
  1045.             for (int i = 0; i < iSize; ++i) {
  1046.                 for (int j = 0; j < iSize; ++j)
  1047.                     aadblR[i][j] = i > j ? DotProduct (aadblQ[i], aadblA[j]) : 0.;
  1048.             }

  1050.             return new org.drip.numerical.linearalgebra.QR (aadblQ, aadblR);
  1051.         } catch (java.lang.Exception e) {
  1052.             e.printStackTrace();
  1053.         }

  1055.         return null;
  1056.     }

  1058.     /**
  1059.      * Retrieve the Triangular Type of the Matrix
  1060.      * 
  1061.      * @param aadblA The Input Matrix
  1062.      * @param dblFloor The Floor Level that means "Zero"
  1063.      * 
  1064.      * @return The Triangular Type
  1065.      */

  1067.     public static final int TriangularType (
  1068.         final double[][] aadblA,
  1069.         final double dblFloor)
  1070.     {
  1071.         if (null == aadblA || !org.drip.numerical.common.NumberUtil.IsValid (dblFloor) || dblFloor < 0.)
  1072.             return NON_TRIANGULAR;

  1074.         int iSize = aadblA.length;
  1075.         boolean bLowerTriangular = true;
  1076.         boolean bUpperTriangular = true;

  1078.         if (1 >= iSize || null == aadblA[0] || iSize != aadblA[0].length) return NON_TRIANGULAR;

  1080.         for (int i = 0; i < iSize; ++i) {
  1081.             for (int j = 0; j < iSize; ++j) {
  1082.                 if (i > j) {
  1083.                     if (java.lang.Math.abs (aadblA[i][j]) > dblFloor) bLowerTriangular = false;
  1084.                 } else if (i < j) {
  1085.                     if (java.lang.Math.abs (aadblA[i][j]) > dblFloor) bUpperTriangular = false;
  1086.                 }
  1087.             }
  1088.         }

  1090.         if (bLowerTriangular && bUpperTriangular) return LOWER_AND_UPPER_TRIANGULAR;

  1092.         if (bLowerTriangular && !bUpperTriangular) return LOWER_TRIANGULAR;

  1094.         if (!bLowerTriangular && bUpperTriangular) return UPPER_TRIANGULAR;

  1096.         return NON_TRIANGULAR;
  1097.     }

  1099.     /**
  1100.      * Compute the Rayleigh Quotient given the Matrix and one of its Eigenvector
  1101.      * 
  1102.      * @param matrix The Given Matrix
  1103.      * @param eigenvector The corresponding Eigenvector
  1104.      * 
  1105.      * @return The Computed Rayleigh Quotient
  1106.      * 
  1107.      * @throws java.lang.Exception Thrown if the Inputs are Invalid
  1108.      */

  1110.     public static final double RayleighQuotient (
  1111.         final double[][] matrix,
  1112.         final double[] eigenvector)
  1113.         throws java.lang.Exception
  1114.     {
  1115.         return org.drip.numerical.linearalgebra.Matrix.DotProduct (
  1116.             eigenvector,
  1117.             org.drip.numerical.linearalgebra.Matrix.Product (
  1118.                 matrix,
  1119.                 eigenvector
  1120.             )
  1121.         );
  1122.     }

  1124.     /**
  1125.      * Scale the Entries of the Input Vector by the Factor
  1126.      * 
  1127.      * @param vector The Input Vector
  1128.      * @param scaleFactor The Scale Factor
  1129.      * 
  1130.      * @return The Scaled Matrix
  1131.      */

  1133.     public static final double[] Scale1D (
  1134.         final double[] vector,
  1135.         final double scaleFactor)
  1136.     {
  1137.         if (null == vector || !org.drip.numerical.common.NumberUtil.IsValid (scaleFactor))
  1138.         {
  1139.             return null;
  1140.         }

  1142.         int rowCount = vector.length;
  1143.         double[] scaledVector = 0 == rowCount ? null : new double[rowCount];

  1145.         for (int rowIndex = 0; rowIndex < rowCount ; ++rowIndex)
  1146.         {
  1147.             scaledVector[rowIndex] = vector[rowIndex] * scaleFactor;
  1148.         }

  1150.         return scaledVector;
  1151.     }

  1153.     /**
  1154.      * Scale the Entries of the Input Matrix by the Factor
  1155.      * 
  1156.      * @param matrix The Input Matrix
  1157.      * @param scaleFactor The Scale Factor
  1158.      * 
  1159.      * @return The Scaled Matrix
  1160.      */

  1162.     public static final double[][] Scale2D (
  1163.         final double[][] matrix,
  1164.         final double scaleFactor)
  1165.     {
  1166.         if (null == matrix || !org.drip.numerical.common.NumberUtil.IsValid (scaleFactor))
  1167.         {
  1168.             return null;
  1169.         }

  1171.         int rowCount = matrix.length;
  1172.         int columnCount = 0 == rowCount || null == matrix[0] ? 0 : matrix[0].length;
  1173.         double[][] scaledMatrix = 0 == columnCount ? null : new double[rowCount][columnCount];

  1175.         for (int rowIndex = 0; rowIndex < rowCount ; ++rowIndex)
  1176.         {
  1177.             for (int columnIndex = 0; columnIndex < columnCount ; ++columnIndex)
  1178.             {
  1179.                 scaledMatrix[rowIndex][columnIndex] = matrix[rowIndex][columnIndex] * scaleFactor;
  1180.             }
  1181.         }

  1183.         return scaledMatrix;
  1184.     }

  1186.     /**
  1187.      * Indicate if the Specified Matrix is Diagonal
  1188.      *
  1189.      * @param matrix The Matrix
  1190.      *
  1191.      * @return TRUE - The Specified Matrix is Diagonal
  1192.      */

  1194.     public static final boolean IsDiagonal (
  1195.         final double[][] matrix)
  1196.     {
  1197.        if (null == matrix)
  1198.        {
  1199.            return false;
  1200.        }

  1202.        int rowCount = matrix.length;
  1203.        int columnCount = 0 == rowCount || null == matrix[0] ? 0 : matrix[0].length;

  1205.        for (int rowIndex = 0;
  1206.             rowIndex < rowCount;
  1207.             ++rowIndex)
  1208.        {
  1209.            for (int columnIndex = 0;
  1210.                 columnIndex < columnCount;
  1211.                 ++columnIndex)
  1212.            {
  1213.                if (rowIndex != columnIndex && 0. != matrix[rowIndex][columnIndex])
  1214.                {
  1215.                    return false;
  1216.                }
  1217.            }
  1218.        }

  1220.        return true;
  1221.     }
  1222. }
Created with JaCoCo 0.7.9.201702052155