C# SV Decomp Example

using System;

using CenterSpace.NMath.Core;
using CenterSpace.NMath.Matrix;

namespace CenterSpace.NMath.Matrix.Examples.CSharp
{
  /// <summary>
  /// A .NET example in C# demonstrating the features of the singular value decomposition (SVD) classes.
  /// </summary>
  class SVDecompExample
  {

    static void Main(string[] args)
    {
      // A general m x n system with random entries.
      RandGenUniform rng = new RandGenUniform(-1, 1);
      rng.Reset(0x124);
      int rows = 6;
      int cols = 3;
      DoubleComplexMatrix A = new DoubleComplexMatrix(rows, cols, rng);

      // Construct a SV decomposition of A.
      DoubleComplexSVDecomp decomp = new DoubleComplexSVDecomp(A);

      // Look at the components of the factorization A = USV'.
      Console.WriteLine();
      Console.WriteLine("U = {0}", decomp.LeftVectors.ToString("F4"));

      Console.WriteLine();
      Console.WriteLine("V = {0}", decomp.RightVectors.ToString("F4"));

      Console.WriteLine();
      Console.WriteLine("s = {0}", decomp.SingularValues);

      // U = 6x3 [ (-0.1773,0.4673)  (-0.2035,0.0897) (0.0996,-0.1514)  
      //           (0.1020,0.0811)   (0.2107,0.1197)  (0.0504,-0.2434) 
      //           (0.0774,-0.5596)  (-0.2614,0.0474) (0.1327,0.3823)  
      //           (0.2475,-0.1521)  (0.0260,0.2871)  (-0.4492,-0.1731) 
      //           (0.1298,-0.5129)  (0.1649,-0.2407) (0.2869,-0.6443)  
      //           (-0.1468,-0.1681) (0.2985,0.7510)  (-0.0955,-0.0581) ]
      //
      // V = 3x3 [ (0.5561,0.0000)  (0.6742,0.0000)   (-0.4859,0.0000) 
      //           (-0.6442,-0.1739) (0.5122,-0.2140) (-0.0265,-0.4959) 
      //           (-0.4539,0.1988)  (0.3357,0.3529)  (-0.0537,0.7172) ]
      //
      // s = [ 2.30473458507626 1.59555281597513 1.05912909184625 ]
      //
      // Note that the singular values, elements on the main diagonal of 
      // the diagonal matrix S, are returned as a vector.

      // The class DoubleComplexSVDecompServer allows more control over the 
      // computation. Suppose that you are only interested in the singular values,
      // not the vectors. You can configure a DoubleComplexSVDecompServer object
      // to compute just the singular values.
      DoubleComplexSVDecompServer decompServer = new DoubleComplexSVDecompServer();
      decompServer.ComputeLeftVectors = false;
      decompServer.ComputeRightVectors = false;
      decomp = decompServer.GetDecomp(A);

      Console.WriteLine();
      Console.WriteLine("Number of left vectors computed: {0}", decomp.NumberLeftVectors); // 0

      Console.WriteLine();
      Console.WriteLine("Number of right vectors computed: {0}", decomp.NumberRightVectors); // 0

      // By default, the "reduced" SVD is computed; that is, if A is m x n, then U
      // is m x n. The "full" SVD is obtained by  adjoining an additional m-n
      // (assuming m > n) orthonormal columns to U making it a m x m unitary matrix. 
      // The singular value decomposition server object can be configured to
      // compute the full SVD as follows:
      decompServer.ComputeLeftVectors = true;
      decompServer.ComputeFull = true;
      decomp = decompServer.GetDecomp(A);

      Console.WriteLine();
      Console.WriteLine("Full SVD, U = {0}", decomp.LeftVectors.ToString("F4"));

      // U = 6x6 
      // [ (-0.1773,0.4673) (-0.2035,0.0897) (0.0996,-0.1514)  (0.1572,-0.1525)  (0.6281,0.4364)  (-0.0526,-0.1796)
      //   (0.1020,0.0811)  (0.2107,0.1197)  (0.0504,-0.2434)  (-0.1824,0.2282)  (-0.0817,0.3425) (-0.6056,0.5352)
      //   (0.0774,-0.5596) (-0.2614,0.0474) (0.1327,0.3823)   (0.0384,0.4186)   (0.4310,0.0384)  (-0.2681,-0.1039)
      //   (0.2475,-0.1521) (0.0260,0.2871)  (-0.4492,-0.1731) (0.7009,-0.0308)  (-0.1625,0.1318) (-0.1724,-0.1873)
      //   (0.1298,-0.5129) (0.1649,-0.2407) (0.2869,-0.6443)  (-0.0732,-0.2579) (0.1492,0.0899)  (0.0388,-0.1840)  
      //   (-0.1468,-0.1681) (0.2985,0.7510) (-0.0955,-0.0581) (-0.3352,0.1166)  (0.0718,0.1567)  (0.3387,-0.1195) ]

      // You can also set a tolerance for the singular values. Singular values 
      // whose value is less than the tolerance are set to zero. The number of
      // singular vectors are adjusted accordingly.
      // 
      // Make A rank deficient.
      A.Col(0)[Slice.All] = A.Col(1); // Two equal columns.
      decompServer.ComputeFull = false;
      decompServer.ComputeLeftVectors = true;
      decompServer.ComputeRightVectors = true;
      decomp = decompServer.GetDecomp(A);

      Console.WriteLine();
      Console.WriteLine("Rank of A = {0}", decomp.Rank); // 3

      // Apparently A has full rank. Let's look at the smallest
      // singular value. Singular values are arranged in descending
      // order, so the smallest value is the last value.
      Console.WriteLine();
      Console.WriteLine("Smallest singular value = {0}", decomp.SingularValue(2)); //  5.38294957152069E-17

      // This singular value is equal to 0, which is within machine precision. Truncating 
      // the SVD will set this value to 0.
      double tolerance = 1e-16;
      decompServer.Tolerance = 1e-16;
      decomp = decompServer.GetDecomp(A);

      // Now look at the rank.
      Console.WriteLine("Rank of A = {0}", decomp.Rank); // 2

      // You can also truncate an existing decomp by calling its
      // truncate method with a specified tolerance.
      tolerance = 1e-12;
      decomp.Truncate(tolerance);

      Console.WriteLine();
      Console.WriteLine("Press Enter Key");
      Console.Read();
    }
  }
}