ooples
diff --git a/‎src/DecompositionMethods/MatrixDecomposition/BidiagonalDecomposition.cs‎
Lines changed: 123 additions & 36 deletions b/‎src/DecompositionMethods/MatrixDecomposition/BidiagonalDecomposition.cs‎
Lines changed: 123 additions & 36 deletions
@@ -1,50 +1,98 @@
 namespace AiDotNet.DecompositionMethods.MatrixDecomposition;
 
-public class BidiagonalDecomposition<T> : IMatrixDecomposition<T>
+/// <summary>
+/// Implements the Bidiagonal Decomposition of a matrix, which factors a matrix into U*B*V^T,
+/// where U and V are orthogonal matrices and B is a bidiagonal matrix.
+/// </summary>
+/// <typeparam name="T">The numeric type used in the matrix (e.g., double, float).</typeparam>
+/// <remarks>
+/// <para>
+/// Bidiagonal decomposition transforms a matrix A into the form A = U*B*V^T, where:
+/// - U is an orthogonal matrix (left singular vectors)
+/// - B is a bidiagonal matrix (non-zero elements only on main diagonal and superdiagonal)
+/// - V^T is the transpose of an orthogonal matrix (right singular vectors)
+/// </para>
+/// <para>
+/// <b>For Beginners:</b> Bidiagonal decomposition is like organizing a complex matrix into a simpler form
+/// that's easier to work with. Think of it like arranging books on a shelf: instead of having them scattered,
+/// you organize them so most spaces are empty (zeros) and important information (non-zero values) is concentrated
+/// in just two lines. This makes many calculations much faster and more efficient, especially when computing
+/// singular values or solving systems of equations.
+/// </para>
+/// <para>
+/// Common applications include:
+/// <list type="bullet">
+/// <item>Computing Singular Value Decomposition (SVD)</item>
+/// <item>Solving least squares problems</item>
+/// <item>Computing eigenvalues of symmetric matrices</item>
+/// <item>Data compression and dimensionality reduction</item>
+/// </list>
+/// </para>
+/// </remarks>
+public class BidiagonalDecomposition<T> : MatrixDecompositionBase<T>
 {
-    private readonly INumericOperations<T> _numOps;
-
-    /// <summary>
-    /// Gets the original matrix that was decomposed.
-    /// </summary>
-    public Matrix<T> A { get; }
-
     /// <summary>
     /// Gets the left orthogonal matrix in the decomposition.
-    /// In simpler terms, this matrix helps transform the original matrix's columns.
     /// </summary>
+    /// <remarks>
+    /// <b>For Beginners:</b> This matrix helps transform the original matrix's columns.
+    /// It represents how the original data is rotated or reflected in the column space.
+    /// </remarks>
     public Matrix<T> U { get; private set; }
 
     /// <summary>
     /// Gets the bidiagonal matrix in the decomposition.
-    /// A bidiagonal matrix is a special matrix where non-zero values appear only on the main diagonal
-    /// and the diagonal immediately above it (called the superdiagonal).
     /// </summary>
+    /// <remarks>
+    /// <b>For Beginners:</b> A bidiagonal matrix is a special matrix where non-zero values appear only
+    /// on the main diagonal and the diagonal immediately above it (called the superdiagonal).
+    /// All other elements are zero, making computations much more efficient.
+    /// </remarks>
     public Matrix<T> B { get; private set; }
 
     /// <summary>
     /// Gets the right orthogonal matrix in the decomposition.
-    /// In simpler terms, this matrix helps transform the original matrix's rows.
     /// </summary>
+    /// <remarks>
+    /// <b>For Beginners:</b> This matrix helps transform the original matrix's rows.
+    /// It represents how the original data is rotated or reflected in the row space.
+    /// </remarks>
     public Matrix<T> V { get; private set; }
 
+    private BidiagonalAlgorithmType _algorithm;
+
     /// <summary>
     /// Creates a new bidiagonal decomposition of the specified matrix.
     /// </summary>
     /// <param name="matrix">The matrix to decompose.</param>
     /// <param name="algorithm">The algorithm to use for decomposition (default is Householder).</param>
     /// <remarks>
+    /// <para>
     /// Bidiagonal decomposition breaks down a matrix A into three simpler matrices: U, B, and V,
     /// where A = U*B*V^T. This makes many matrix operations easier to perform.
+    /// </para>
+    /// <para>
+    /// <b>For Beginners:</b> This constructor takes your original matrix and breaks it down into simpler
+    /// components using the specified algorithm. The Householder algorithm (default) is generally the most
+    /// stable and efficient method for most matrices.
+    /// </para>
     /// </remarks>
     public BidiagonalDecomposition(Matrix<T> matrix, BidiagonalAlgorithmType algorithm = BidiagonalAlgorithmType.Householder)
+        : base(matrix)
     {
-        A = matrix;
-        _numOps = MathHelper.GetNumericOperations<T>();
+        _algorithm = algorithm;
         U = new Matrix<T>(matrix.Rows, matrix.Rows);
         B = new Matrix<T>(matrix.Columns, matrix.Columns);
         V = new Matrix<T>(matrix.Columns, matrix.Columns);
-        Decompose(algorithm);
+        Decompose();
+    }
+
+    /// <summary>
+    /// Performs the bidiagonal decomposition.
+    /// </summary>
+    protected override void Decompose()
+    {
+        ComputeDecomposition(_algorithm);
     }
 
     /// <summary>
@@ -57,6 +105,17 @@ public BidiagonalDecomposition(Matrix<T> matrix, BidiagonalAlgorithmType algorit
     /// Householder is generally the most stable and commonly used method.
     /// </remarks>
     public void Decompose(BidiagonalAlgorithmType algorithm = BidiagonalAlgorithmType.Householder)
+    {
+        _algorithm = algorithm;
+        ComputeDecomposition(algorithm);
+    }
+
+    /// <summary>
+    /// Performs the actual decomposition computation using the specified algorithm.
+    /// </summary>
+    /// <param name="algorithm">The algorithm to use for decomposition.</param>
+    /// <exception cref="ArgumentException">Thrown when an unsupported algorithm is specified.</exception>
+    private void ComputeDecomposition(BidiagonalAlgorithmType algorithm)
     {
         switch (algorithm)
         {
@@ -96,7 +155,7 @@ private void DecomposeHouseholder()
             Vector<T> v = HouseholderVector(x);
 
             // Apply Householder reflection to B
-            Matrix<T> P = Matrix<T>.CreateIdentity(m - k).Subtract(v.OuterProduct(v).Multiply(_numOps.FromDouble(2)));
+            Matrix<T> P = Matrix<T>.CreateIdentity(m - k).Subtract(v.OuterProduct(v).Multiply(NumOps.FromDouble(2)));
             Matrix<T> subB = B.GetSubMatrix(k, k, m - k, n - k);
             B.SetSubMatrix(k, k, P.Multiply(subB));
 
@@ -111,7 +170,7 @@ private void DecomposeHouseholder()
                 v = HouseholderVector(x);
 
                 // Apply Householder reflection to B
-                P = Matrix<T>.CreateIdentity(n - k - 1).Subtract(v.OuterProduct(v).Multiply(_numOps.FromDouble(2)));
+                P = Matrix<T>.CreateIdentity(n - k - 1).Subtract(v.OuterProduct(v).Multiply(NumOps.FromDouble(2)));
                 subB = B.GetSubMatrix(k, k + 1, m - k, n - k - 1);
                 B.SetSubMatrix(k, k + 1, subB.Multiply(P));
 
@@ -177,7 +236,7 @@ private void DecomposeLanczos()
         Random rand = new();
         for (int i = 0; i < n; i++)
         {
-            v[i] = _numOps.FromDouble(rand.NextDouble());
+            v[i] = NumOps.FromDouble(rand.NextDouble());
         }
         v = v.Divide(v.Norm());
 
@@ -206,10 +265,17 @@ private void DecomposeLanczos()
     /// <returns>The solution vector x.</returns>
     /// <exception cref="ArgumentException">Thrown when the length of vector b doesn't match the number of rows in matrix A.</exception>
     /// <remarks>
+    /// <para>
     /// This method uses the decomposition to efficiently solve the system without
     /// directly inverting the matrix, which is more numerically stable.
+    /// </para>
+    /// <para>
+    /// <b>For Beginners:</b> This finds the values of x that satisfy Ax = b. Using the bidiagonal
+    /// decomposition makes this much faster than directly solving the original equation, especially
+    /// for large matrices. It's like solving a puzzle by first organizing the pieces into groups.
+    /// </para>
     /// </remarks>
-    public Vector<T> Solve(Vector<T> b)
+    public override Vector<T> Solve(Vector<T> b)
     {
         if (b.Length != A.Rows)
             throw new ArgumentException("Vector b must have the same length as the number of rows in matrix A.");
@@ -225,18 +291,25 @@ public Vector<T> Solve(Vector<T> b)
     /// </summary>
     /// <returns>The inverse matrix of A.</returns>
     /// <remarks>
+    /// <para>
     /// Matrix inversion is computationally expensive and can be numerically unstable.
     /// When possible, use the Solve method instead of explicitly computing the inverse.
+    /// </para>
+    /// <para>
+    /// <b>For Beginners:</b> The inverse of a matrix is like the reciprocal of a number.
+    /// Just as 5 * (1/5) = 1, a matrix multiplied by its inverse gives the identity matrix.
+    /// However, computing the inverse is expensive, so it's usually better to use Solve() instead.
+    /// </para>
     /// </remarks>
-    public Matrix<T> Invert()
+    public override Matrix<T> Invert()
     {
         int n = A.Columns;
         Matrix<T> inverse = new Matrix<T>(n, n);
 
         for (int i = 0; i < n; i++)
         {
             Vector<T> ei = new Vector<T>(n);
-            ei[i] = _numOps.One;
+            ei[i] = NumOps.One;
             inverse.SetColumn(i, Solve(ei));
         }
 
@@ -249,19 +322,26 @@ public Matrix<T> Invert()
     /// <param name="x">The input vector.</param>
     /// <returns>A Householder vector that can be used to create a reflection matrix.</returns>
     /// <remarks>
+    /// <para>
     /// A Householder reflection is a transformation that reflects a vector across a plane.
     /// It's used to zero out specific elements in a matrix during decomposition.
+    /// </para>
+    /// <para>
+    /// <b>For Beginners:</b> Think of a Householder reflection like holding up a mirror to a vector.
+    /// The reflection helps us systematically create zeros in the matrix, simplifying its structure
+    /// step by step during the decomposition process.
+    /// </para>
     /// </remarks>
     private Vector<T> HouseholderVector(Vector<T> x)
     {
         T norm = x.Norm();
         Vector<T> v = x.Clone();
-        v[0] = _numOps.Add(v[0], _numOps.Multiply(_numOps.SignOrZero(x[0]), norm));
+        v[0] = NumOps.Add(v[0], NumOps.Multiply(NumOps.SignOrZero(x[0]), norm));
 
         return v.Divide(v.Norm());
     }
 
-        /// <summary>
+    /// <summary>
     /// Applies a Givens rotation to the specified matrices.
     /// </summary>
     /// <param name="M">The matrix to which the rotation is applied.</param>
@@ -272,33 +352,40 @@ private Vector<T> HouseholderVector(Vector<T> x)
     /// <param name="l">The second column/row index for the rotation.</param>
     /// <param name="isLeft">If true, applies a left rotation (row operation); otherwise, applies a right rotation (column operation).</param>
     /// <remarks>
+    /// <para>
     /// A Givens rotation is a simple rotation in a plane spanned by two coordinate axes.
     /// It's used to selectively zero out specific elements in a matrix during decomposition.
+    /// </para>
+    /// <para>
+    /// <b>For Beginners:</b> A Givens rotation is like turning a dial to make specific values become zero.
+    /// It's precise and only affects two rows or columns at a time, making it useful for sparse matrices
+    /// or when you need to zero out specific elements without disturbing others.
+    /// </para>
     /// </remarks>
     private void GivensRotation(Matrix<T> M, Matrix<T> Q, int i, int k, int j, int l, bool isLeft)
     {
         T a = M[i, j];
         T b = M[k, l];
-        T r = _numOps.Sqrt(_numOps.Add(_numOps.Multiply(a, a), _numOps.Multiply(b, b)));
-        T c = _numOps.Divide(a, r);
-        T s = _numOps.Divide(b, r);
+        T r = NumOps.Sqrt(NumOps.Add(NumOps.Multiply(a, a), NumOps.Multiply(b, b)));
+        T c = NumOps.Divide(a, r);
+        T s = NumOps.Divide(b, r);
 
         if (isLeft)
         {
             for (int j2 = 0; j2 < M.Columns; j2++)
             {
                 T temp1 = M[i, j2];
                 T temp2 = M[k, j2];
-                M[i, j2] = _numOps.Add(_numOps.Multiply(c, temp1), _numOps.Multiply(s, temp2));
-                M[k, j2] = _numOps.Subtract(_numOps.Multiply(_numOps.Negate(s), temp1), _numOps.Multiply(c, temp2));
+                M[i, j2] = NumOps.Add(NumOps.Multiply(c, temp1), NumOps.Multiply(s, temp2));
+                M[k, j2] = NumOps.Subtract(NumOps.Multiply(NumOps.Negate(s), temp1), NumOps.Multiply(c, temp2));
             }
 
             for (int i2 = 0; i2 < Q.Rows; i2++)
             {
                 T temp1 = Q[i2, i];
                 T temp2 = Q[i2, k];
-                Q[i2, i] = _numOps.Add(_numOps.Multiply(c, temp1), _numOps.Multiply(s, temp2));
-                Q[i2, k] = _numOps.Subtract(_numOps.Multiply(_numOps.Negate(s), temp1), _numOps.Multiply(c, temp2));
+                Q[i2, i] = NumOps.Add(NumOps.Multiply(c, temp1), NumOps.Multiply(s, temp2));
+                Q[i2, k] = NumOps.Subtract(NumOps.Multiply(NumOps.Negate(s), temp1), NumOps.Multiply(c, temp2));
             }
         }
         else
@@ -307,16 +394,16 @@ private void GivensRotation(Matrix<T> M, Matrix<T> Q, int i, int k, int j, int l
             {
                 T temp1 = M[i2, j];
                 T temp2 = M[i2, l];
-                M[i2, j] = _numOps.Add(_numOps.Multiply(c, temp1), _numOps.Multiply(s, temp2));
-                M[i2, l] = _numOps.Subtract(_numOps.Multiply(_numOps.Negate(s), temp1), _numOps.Multiply(c, temp2));
+                M[i2, j] = NumOps.Add(NumOps.Multiply(c, temp1), NumOps.Multiply(s, temp2));
+                M[i2, l] = NumOps.Subtract(NumOps.Multiply(NumOps.Negate(s), temp1), NumOps.Multiply(c, temp2));
             }
 
             for (int j2 = 0; j2 < Q.Columns; j2++)
             {
                 T temp1 = Q[j, j2];
                 T temp2 = Q[l, j2];
-                Q[j, j2] = _numOps.Add(_numOps.Multiply(c, temp1), _numOps.Multiply(s, temp2));
-                Q[l, j2] = _numOps.Subtract(_numOps.Multiply(_numOps.Negate(s), temp1), _numOps.Multiply(c, temp2));
+                Q[j, j2] = NumOps.Add(NumOps.Multiply(c, temp1), NumOps.Multiply(s, temp2));
+                Q[l, j2] = NumOps.Subtract(NumOps.Multiply(NumOps.Negate(s), temp1), NumOps.Multiply(c, temp2));
             }
         }
     }
@@ -340,8 +427,8 @@ private Vector<T> SolveBidiagonal(Vector<T> y)
         {
             T sum = y[i];
             if (i < n - 1)
-                sum = _numOps.Subtract(sum, _numOps.Multiply(B[i, i + 1], x[i + 1]));
-            x[i] = _numOps.Divide(sum, B[i, i]);
+                sum = NumOps.Subtract(sum, NumOps.Multiply(B[i, i + 1], x[i + 1]));
+            x[i] = NumOps.Divide(sum, B[i, i]);
         }
 
         return x;