diff --git a/Docs/Internal/NumericalAnalysis/NumericalAnalysis_v6.41.pdf b/Docs/Internal/NumericalAnalysis/NumericalAnalysis_v6.42.pdf similarity index 65% rename from Docs/Internal/NumericalAnalysis/NumericalAnalysis_v6.41.pdf rename to Docs/Internal/NumericalAnalysis/NumericalAnalysis_v6.42.pdf index 77393d6ec10..f4939e1182c 100644 Binary files a/Docs/Internal/NumericalAnalysis/NumericalAnalysis_v6.41.pdf and b/Docs/Internal/NumericalAnalysis/NumericalAnalysis_v6.42.pdf differ diff --git a/Docs/Internal/NumericalAnalysis/Shruti/Source/NumericalAnalysis_v6.41.docx b/Docs/Internal/NumericalAnalysis/Shruti/Source/NumericalAnalysis_v6.41.docx index c28cb52e792..e786c96f078 100644 Binary files a/Docs/Internal/NumericalAnalysis/Shruti/Source/NumericalAnalysis_v6.41.docx and b/Docs/Internal/NumericalAnalysis/Shruti/Source/NumericalAnalysis_v6.41.docx differ diff --git a/Docs/Internal/NumericalAnalysis/Shruti/Source/NumericalAnalysis_v6.42.docx b/Docs/Internal/NumericalAnalysis/Shruti/Source/NumericalAnalysis_v6.42.docx new file mode 100644 index 00000000000..7b48abcab3a Binary files /dev/null and b/Docs/Internal/NumericalAnalysis/Shruti/Source/NumericalAnalysis_v6.42.docx differ diff --git a/IdeaDRIP/NumericalAnalysis/NA_v0.05 b/IdeaDRIP/NumericalAnalysis/NA_v0.05 new file mode 100644 index 00000000000..3d7bda4f743 --- /dev/null +++ b/IdeaDRIP/NumericalAnalysis/NA_v0.05 @@ -0,0 +1,153 @@ + + -------------------------- + #1 - Successive Over-Relaxation + -------------------------- + -------------------------- + 1.1) Successive Over-Relaxation Method for A.x = b; A - n x n square matrix, x - unknown vector, b - RHS vector + 1.2) Decompose A into Diagonal, Lower, and upper Triangular Matrices; A = D + L + U + 1.3) SOR Scheme uses omega input + 1.4) Forward subsitution scheme to iteratively determine the x_i's + 1.5) SOR Scheme Linear System Convergence: Inputs A and D, Jacobi Iteration Matrix Spectral Radius, omega + - Construct Jacobi Iteration Matrix: C_Jacobian = I - (Inverse D) A + - Convergence Verification #1: Ensure that Jacobi Iteration Matrix Spectral Radius is < 1 + - Convergence Verification #2: Ensure omega between 0 and 2 + - Optimal Relaxation Parameter Expression in terms of Jacobi Iteration Matrix Spectral Radius + - Omega Based Convergence Rate Expression + - Gauss-Seidel omega = 1; corresponding Convergence Rate + - Optimal Omega Convergence Rate + 1.6) Generic Iterative Solver Method: + - Inputs: Iterator Function(x) and omega + - Unrelaxed Iterated variable: x_n+1 = f(x_n) + - SOR Based Iterated variable: x_n+1 = (1-omega).x_n + omega.f(x_n) + - SOR Based Iterated variable for Unknown Vector x: x_n+1 = (1-omega).x_n + omega.(L_star inverse)(b - U.x_n) + -------------------------- + + -------------------------- + #2 - Successive Over-Relaxation + -------------------------- + -------------------------- + 2.1) SSOR Algorithm - Inputs; A, omega, and gamma + - Decompose into D and L + - Pre-conditioner Matrix: Expression from SSOR + - Finally SSOR Iteration Formula + -------------------------- + + ---------------------------- + #7 - Tridiagonal matrix algorithm + ---------------------------- + ---------------------------- + 7.1) Is Tridiagonal Check + 7.2) Core Algorithm: + - C Prime's and D Prime's Calculation + - Back Substitution for the Result + - Modified better Book-keeping algorithm + 7.3) Sherman-Morrison Algorithm: + - Choice of gamma + - Construct Tridiagonal B from A and gamma + - u Column from gamma and c_n + - v Column from a_1 and gamma + - Solve for y from By=d + - Solve for q from Bq=u + - Use Sherman Morrison Formula to extract x + 7.4) Alternate Boundary Condition Algorithm: + - Solve Au=d for u + - Solve Av={-a_2, 0, ..., -c_n} for v + - Full solution is x_i = u_i + x_1 * v_i + - x_1 if computed using formula + ---------------------------- + + ---------------------------- + #8 - Triangular Matrix + ---------------------------- + ---------------------------- + 8.1) Description: + - Lower/Left Triangular Verification + - Upper/Right Triangular Verification + - Diagonal Matrix Verification + - Upper/Lower Trapezoidal Verification + 8.2) Forward/Back Substitution: + - Inputs => L and b + - Forward Substitution + - Inputs => U and b + - Back Substitution + 8.3) Properties: + - Is Matrix Normal, i.e., A times A transpose = A transpose times A + - Characteristic Polynomial + - Determinant/Permanent + 8.4) Special Forms: + - Upper/Lower Unitriangular Matrix Verification + - Upper/Lower Strictly Matrix Verification + - Upper/Lower Atomic Matrix Verification + ---------------------------- + + ---------------------------- + #9 - Sylvester Equation + ---------------------------- + ---------------------------- + 9.1) Matrix Form: + - Inputs: A, B, and C + - Size Constraints Verification + 9.2) Solution Criteria: + - Co-joint EigenSpectrum between A and B + 9.3) Numerical Solution: + - Decomposition of A/B using Schur Decomposition into Triangular Form + - Forward/Back Substitution + ---------------------------- + + ---------------------------- + #10 - Bartels-Stewart Algorithm + ---------------------------- + ---------------------------- + 10.1) Matrix Form: + - Inputs: A, B, and C + - Size Constraints Verification + 10.2) Schur Decompositions: + - R = U^T A U - emits U and R + - S = V^T B^T V - emits V and S + - F = U^T C V + - Y = U^T X V + - Solution to R.Y - Y.S^T = F + - Finally X = U.Y.V^T + 10.3) Computational Costs: + - Flops cost for Schur decomposition + - Flops cost overall + 10.4) Hessenberg-Schur Decompositions: + - R = U^T A U becomes H = Q^T A Q - thus emits Q and H (Upper Hessenberg) + - Computational Costs + ---------------------------- + + -------------------------- + #12 - Crank–Nicolson method + -------------------------- + -------------------------- + 12.1) von Neumann Stability Validator - Inputs; time-step, diffusivity, space step + - 1D => time step * diffusivity / space step ^2 < 0.5 + - nD => time step * diffusivity / (space step hypotenuse) ^ 2 < 0.5 + 12.2) Set up: + - Input: Spatial variable x + - Input: Time variable t + - Inputs: State variable u, du/dx, d2u/dx2 - all at time t + - Second Order, 1D => du/dt = F (u, x, t, du/dx, d2u/dx2) + 12.3) Finite Difference Evolution Scheme: + - Time Step delta_t, space step delta_x + - Forward Difference: F calculated at x + - Backward Difference: F calculated at x + delta_x + - Crank Nicolson: Average of Forward/Backward + 12.4) 1D Diffusion: + - Inputs: 1D von Neumann Stability Validator, Number of time/space steps + - Time Index n, Space Index i + - Explicit Tridiagonal form for the discretization - State concentration at n+1 given state concentration at n + - Non-linear diffusion coefficient: + - Linearization across x_i and x_i+1 + - Quasi-explicit forms accommodation + 12.5) 2D Diffusion: + - Inputs: 2D von Neumann Stability Validator, Number of time/space steps + - Time Index n, Space Index i, j + - Explicit Tridiagonal form for the discretization - State concentration at n+1 given state concentration at n + - Explicit Solution using the Alternative Difference Implicit Scheme + 12.6) Extension to Iterative Solver Schemes above: + - Input: State Space "velocity" dF/du + - Input: State "Step Size" delta_u + - Fixed Point Iterative Location Scheme + - Relaxation Scheme based Robustness => Input: Relaxation Parameter + -------------------------- diff --git a/NumericalAnalysisLibrary.md b/NumericalAnalysisLibrary.md index b3fcebdce38..92ed2249a09 100644 --- a/NumericalAnalysisLibrary.md +++ b/NumericalAnalysisLibrary.md @@ -349,12 +349,19 @@ Numerical Analysis Library contains the supporting Functionality for Numerical M * Crank–Nicolson for nonlinear problems * Application in financial mathematics * References - * Sylvester equation + * Sylvester Equation * Existence and uniqueness of the solutions * Roth's removal rule * Numerical solutions * References - * Crank–Nicolson Method + * Bartels–Stewart Algorithm + * Introduction + * The Algorithm + * Computational Cost + * Simplifications and Special Cases + * The Hessenberg–Schur algorithm + * References + * Triangular Matrix * Description * Forward and Back Substitution * Forward Substitution diff --git a/ReleaseNotes/02_14_2024.txt b/ReleaseNotes/02_14_2024.txt new file mode 100644 index 00000000000..32bc327b41d --- /dev/null +++ b/ReleaseNotes/02_14_2024.txt @@ -0,0 +1,25 @@ + +Features: + +Bug Fixes/Re-organization: + +Samples: + +IdeaDRIP: + + - Bartels-Stewart Algorithm The - Simplifications and Special Cases (1) + - Bartels-Stewart Algorithm The Hessenberg–Schur Algorithm (2-8) + - Bartels-Stewart Algorithm Software and Implementation (9-10) + - Bartels-Stewart Algorithm Alternate Approaches (11-15) + - Alternating-direction Implicit Method Introduction (16-21) + - Alternating-direction Implicit Method ADI for Matrix Equations - The (22-29) + - Alternating-direction Implicit Method ADI for Matrix Equations - When to use (30-47) + - Alternating-direction Implicit Method ADI for Matrix Equations - Shift-parameter Selection and the Error Equation (48-57) + - Alternating-direction Implicit Method ADI for Matrix Equations - Shift-parameter Selection and the Error Equation - Near-optimal Shift Parameters (58-67) + - Alternating-direction Implicit Method ADI for Matrix Equations - Shift-parameter Selection and the Error Equation - Heuristic Strategies (68-73) + - Alternating-direction Implicit Method ADI for Matrix Equations - Factored (74-79) + - Alternating-direction Implicit Method ADI for Parabolic Equations (80-85) + - Alternating-direction Implicit Method ADI for Parabolic Equations - 2D Diffusion Equations (86-107) + - Alternating-direction Implicit Method ADI for Parabolic Equations - Generalizations (108-109) + - Alternating-direction Implicit Method Fundamental ADI (FADI) - Simplification of to (110-115) + - Alternating-direction Implicit Method Fundamental ADI (FADI) - Relations to other methods (116-120) diff --git a/ReleaseNotes/02_15_2024.txt b/ReleaseNotes/02_15_2024.txt new file mode 100644 index 00000000000..67e5f8bdd0d --- /dev/null +++ b/ReleaseNotes/02_15_2024.txt @@ -0,0 +1,45 @@ + +Features: + + - Bartels-Stewart Algorithm Diagnostics #1 (63, 64) + - Bartels-Stewart Algorithm Diagnostics #2 (65, 66) + - Sylvester Equation Solution Annotated Shell (67, 68, 69) + + +Bug Fixes/Re-organization: + + - Numerical Analyzer Bartels-Stewart Algorithm (62) + + +Samples: + + - QR Decomposition Matrix Util (63, 64) + - Graham Schmidt Matrix Orthonormalization Process (70, 71, 72) + - Bartels-Stewart Algorithm Run #1 (73, 74, 75) + - Bartels-Stewart Algorithm Run #2 (76, 77, 78) + - Bartels-Stewart Algorithm Run #3 (79, 80) + - Bartels-Stewart Algorithm Run #4 (81, 82) + - Bartels-Stewart Algorithm Run #5 (83, 84) + - Bartels-Stewart Algorithm Run #6 (85, 86, 87) + - Bartels-Stewart Algorithm Run #7 (88, 89, 90) + - Bartels-Stewart Algorithm Run #8 (91, 92, 93) + - Bartels-Stewart Algorithm Run #9 (94, 95, 96) + - Bartels-Stewart Algorithm Run #10 (97, 98, 99) + - Bartels-Stewart Algorithm Run #11 (100, 101, 102) + - Bartels-Stewart Algorithm Run #12 (103, 104, 105) + - Bartels-Stewart Algorithm Run #13 (106, 107, 108) + - Bartels-Stewart Algorithm Run #14 (109, 110, 111) + - Bartels-Stewart Algorithm Run #15 (112, 113, 114) + - Bartels-Stewart Algorithm Run #16 (115, 116, 117) + - Bartels-Stewart Algorithm Run #17 (118, 119, 120) + + +IdeaDRIP: + + - Alternating-direction Implicit Method Fundamental ADI (FADI) - Relations to other methods (1) + - Bartels-Stewart Algorithm #1 (2-7) + - The Bartels-Stewart Algorithm #1 (8-27) + - Bartels-Stewart Algorithm Hessenberg-Schur #1 (28-31) + - Bartels-Stewart Algorithm #2 (32-37) + - The Bartels-Stewart Algorithm #2 (38-57) + - Bartels-Stewart Algorithm Hessenberg-Schur #2 (58-61) diff --git a/src/main/java/org/drip/numerical/linearsolver/BartelsStewartScheme.java b/src/main/java/org/drip/numerical/linearsolver/BartelsStewartScheme.java index 71b671e97f6..dc708885b47 100644 --- a/src/main/java/org/drip/numerical/linearsolver/BartelsStewartScheme.java +++ b/src/main/java/org/drip/numerical/linearsolver/BartelsStewartScheme.java @@ -126,6 +126,7 @@ public class BartelsStewartScheme { private double[][] _rhsMatrix = null; + private boolean _diagnosticsOn = false; private SylvesterEquation _sylvesterEquation = null; /** @@ -133,13 +134,15 @@ public class BartelsStewartScheme * * @param sylvesterEquation Sylvester Equation Instance * @param rhsMatrix "RHS" Matrix + * @param diagnosticsOn Diagnostics-on Flag * * @throws Exception Thrown if the Inputs are Invalid */ public BartelsStewartScheme ( final SylvesterEquation sylvesterEquation, - final double[][] rhsMatrix) + final double[][] rhsMatrix, + final boolean diagnosticsOn) throws Exception { if (null == (_sylvesterEquation = sylvesterEquation) || @@ -168,6 +171,17 @@ public SylvesterEquation sylvesterEquation() * @return "RHS" Matrix */ + public boolean diagnosticsOn() + { + return _diagnosticsOn; + } + + /** + * Retrieve the Diagnostics On Flag + * + * @return Diagnostics On Flag + */ + public double[][] rhsMatrix() { return _rhsMatrix; @@ -179,20 +193,94 @@ public void solve() System.out.println(); - for (int i = 0; i < qrA.q().length; ++i) { + double[][] u = qrA.q(); + + for (int i = 0; i < u.length; ++i) { + System.out.println ( + "\t| Matrix U => [" + NumberUtil.ArrayRow (u[i], 2, 4, false) + " ]||" + ); + } + + System.out.println(); + + double[][] r = qrA.r(); + + for (int i = 0; i < r.length; ++i) { + System.out.println ( + "\t| Matrix R => [" + NumberUtil.ArrayRow (r[i], 2, 4, false) + " ]||" + ); + } + + System.out.println(); + + QR qrBTranspose = MatrixUtil.QRDecomposition ( + _sylvesterEquation.squareMatrixB().transpose().r2Array() + ); + + double[][] v = qrBTranspose.q(); + + for (int i = 0; i < v.length; ++i) { System.out.println ( - "\t| Matrix Q => [" + NumberUtil.ArrayRow (qrA.q()[i], 2, 4, false) + " ]||" + "\t| Matrix V => [" + NumberUtil.ArrayRow (v[i], 2, 4, false) + " ]||" ); } System.out.println(); - for (int i = 0; i < qrA.r().length; ++i) { + double[][] s = qrBTranspose.r(); + + for (int i = 0; i < s.length; ++i) { System.out.println ( - "\t| Matrix R => [" + NumberUtil.ArrayRow (qrA.r()[i], 2, 4, false) + " ]||" + "\t| Matrix S => [" + NumberUtil.ArrayRow (s[i], 2, 4, false) + " ]||" ); } + double[][] f = MatrixUtil.Product (MatrixUtil.Transpose (u), _rhsMatrix); + + f = MatrixUtil.Product (f, v); + + System.out.println(); + + for (int i = 0; i < f.length; ++i) { + System.out.println ( + "\t| Matrix F => [" + NumberUtil.ArrayRow (f[i], 2, 4, false) + " ]||" + ); + } + + double[][] y = new double[f.length][f.length]; + + for (int i = 0; i < f.length; ++i) { + for (int j = 0; j < f.length; ++j) { + y[i][j] = 0.; + } + } + + for (int i = f.length - 1; i >= 0 ; --i) { + for (int j = f.length - 1; j >= 0; --j) { + double coefficientIJ = 0.; + double ryCumulativeSum = 0.; + double syCumulativeSum = 0.; + + for (int k = i; k < f.length; ++k) { + if (k == i) { + coefficientIJ += r[i][k]; + } else { + ryCumulativeSum += r[i][k] * y[k][j]; + } + } + + for (int k = j; k < f.length; ++k) { + if (k == j) { + coefficientIJ += r[i][k]; + } else { + syCumulativeSum += r[i][k] * y[k][j]; + } + } + + y[i][j] = (f[i][j] - ryCumulativeSum - syCumulativeSum) / coefficientIJ; + } + } + System.out.println(); } @@ -219,7 +307,8 @@ public static final void main ( SquareMatrix.Standard (matrixA), SquareMatrix.Standard (matrixB) ), - matrixC + matrixC, + true ); bartelsStewartScheme.solve(); diff --git a/src/main/java/org/drip/numerical/linearsolver/SylvesterEquationSolution.java b/src/main/java/org/drip/numerical/linearsolver/SylvesterEquationSolution.java new file mode 100644 index 00000000000..0654104e8ae --- /dev/null +++ b/src/main/java/org/drip/numerical/linearsolver/SylvesterEquationSolution.java @@ -0,0 +1,122 @@ + +package org.drip.numerical.linearsolver; + +/* + * -*- mode: java; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- + */ + +/*! + * Copyright (C) 2025 Lakshmi Krishnamurthy + * + * This file is part of DROP, an open-source library targeting analytics/risk, transaction cost analytics, + * asset liability management analytics, capital, exposure, and margin analytics, valuation adjustment + * analytics, and portfolio construction analytics within and across fixed income, credit, commodity, + * equity, FX, and structured products. It also includes auxiliary libraries for algorithm support, + * numerical analysis, numerical optimization, spline builder, model validation, statistical learning, + * graph builder/navigator, and computational support. + * + * https://lakshmidrip.github.io/DROP/ + * + * DROP is composed of three modules: + * + * - DROP Product Core - https://lakshmidrip.github.io/DROP-Product-Core/ + * - DROP Portfolio Core - https://lakshmidrip.github.io/DROP-Portfolio-Core/ + * - DROP Computational Core - https://lakshmidrip.github.io/DROP-Computational-Core/ + * + * DROP Product Core implements libraries for the following: + * - Fixed Income Analytics + * - Loan Analytics + * - Transaction Cost Analytics + * + * DROP Portfolio Core implements libraries for the following: + * - Asset Allocation Analytics + * - Asset Liability Management Analytics + * - Capital Estimation Analytics + * - Exposure Analytics + * - Margin Analytics + * - XVA Analytics + * + * DROP Computational Core implements libraries for the following: + * - Algorithm Support + * - Computation Support + * - Function Analysis + * - Graph Algorithm + * - Model Validation + * - Numerical Analysis + * - Numerical Optimizer + * - Spline Builder + * - Statistical Learning + * + * Documentation for DROP is Spread Over: + * + * - Main => https://lakshmidrip.github.io/DROP/ + * - Wiki => https://github.com/lakshmiDRIP/DROP/wiki + * - GitHub => https://github.com/lakshmiDRIP/DROP + * - Repo Layout Taxonomy => https://github.com/lakshmiDRIP/DROP/blob/master/Taxonomy.md + * - Javadoc => https://lakshmidrip.github.io/DROP/Javadoc/index.html + * - Technical Specifications => https://github.com/lakshmiDRIP/DROP/tree/master/Docs/Internal + * - Release Versions => https://lakshmidrip.github.io/DROP/version.html + * - Community Credits => https://lakshmidrip.github.io/DROP/credits.html + * - Issues Catalog => https://github.com/lakshmiDRIP/DROP/issues + * + * 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 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * 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. + */ + +/** + * SylvesterEquationSolution holds the Solution to the Sylvester Equation, which is defined by: + * + * A.X + X.B = C + * + * Here A, B, and C are the Sylvester Equation components, X is the unknown whose solution is to sought. The + * References are: + * + *

+ * + * + *

+ * + *

+ * + * @author Lakshmi Krishnamurthy + */ + +public class SylvesterEquationSolution +{ + +} diff --git a/src/main/java/org/drip/sample/matrix/GrahamSchmidtProcess.java b/src/main/java/org/drip/sample/matrix/GrahamSchmidtProcess.java index 36e61520f76..7491d51bf9d 100644 --- a/src/main/java/org/drip/sample/matrix/GrahamSchmidtProcess.java +++ b/src/main/java/org/drip/sample/matrix/GrahamSchmidtProcess.java @@ -187,11 +187,6 @@ public static final void main ( MatrixUtil.Product (MatrixUtil.Transpose (aadblUOrthonormal), uvTranspose) ); - /* NumberUtil.PrintMatrix ( - "R CHECK #2", - MatrixUtil.Product (aadblV, MatrixUtil.Transpose (aadblUOrthonormal)) - ); */ - EnvManager.TerminateEnv(); } }