How to resolve the algorithm Cholesky decomposition step by step in the C# programming language
How to resolve the algorithm Cholesky decomposition step by step in the C# programming language
Table of Contents
Problem Statement
Every symmetric, positive definite matrix A can be decomposed into a product of a unique lower triangular matrix L and its transpose:
L
{\displaystyle L}
is called the Cholesky factor of
A
{\displaystyle A}
, and can be interpreted as a generalized square root of
A
{\displaystyle A}
, as described in Cholesky decomposition. In a 3x3 example, we have to solve the following system of equations: We can see that for the diagonal elements (
l
k k
{\displaystyle l_{kk}}
) of
L
{\displaystyle L}
there is a calculation pattern: or in general: For the elements below the diagonal (
l
i k
{\displaystyle l_{ik}}
, where
i
k
{\displaystyle i>k}
) there is also a calculation pattern: which can also be expressed in a general formula: Task description The task is to implement a routine which will return a lower Cholesky factor
L
{\displaystyle L}
for every given symmetric, positive definite nxn matrix
A
{\displaystyle A}
. You should then test it on the following two examples and include your output. Example 1: Example 2:
Let's start with the solution:
Step by Step solution about How to resolve the algorithm Cholesky decomposition step by step in the C# programming language
The code starts by importing the necessary libraries, then defines a method called Print that prints the matrix.
Then defines a function called Cholesky that takes as input a matrix and returns the lower Cholesky Factor of that matrix.
The function first checks if the matrix is square, symmetric, and positive definite, and if not, it returns an error.
Then, it creates a new matrix, ret, which will store the lower Cholesky Factor.
The function then iterates over the rows and columns of the input matrix, and for each element, it calculates the corresponding element of the lower Cholesky Factor.
The function returns the lower Cholesky Factor, which is a lower triangular matrix with positive diagonal elements.
The Main function then defines two test matrices, test1 and test2, and calls the Cholesky function on each of them, printing the results.
Source code in the csharp programming language
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace Cholesky
{
class Program
{
/// <summary>
/// This is example is written in C#, and compiles with .NET Framework 4.0
/// </summary>
/// <param name="args"></param>
static void Main(string[] args)
{
double[,] test1 = new double[,]
{
{25, 15, -5},
{15, 18, 0},
{-5, 0, 11},
};
double[,] test2 = new double[,]
{
{18, 22, 54, 42},
{22, 70, 86, 62},
{54, 86, 174, 134},
{42, 62, 134, 106},
};
double[,] chol1 = Cholesky(test1);
double[,] chol2 = Cholesky(test2);
Console.WriteLine("Test 1: ");
Print(test1);
Console.WriteLine("");
Console.WriteLine("Lower Cholesky 1: ");
Print(chol1);
Console.WriteLine("");
Console.WriteLine("Test 2: ");
Print(test2);
Console.WriteLine("");
Console.WriteLine("Lower Cholesky 2: ");
Print(chol2);
}
public static void Print(double[,] a)
{
int n = (int)Math.Sqrt(a.Length);
StringBuilder sb = new StringBuilder();
for (int r = 0; r < n; r++)
{
string s = "";
for (int c = 0; c < n; c++)
{
s += a[r, c].ToString("f5").PadLeft(9) + ",";
}
sb.AppendLine(s);
}
Console.WriteLine(sb.ToString());
}
/// <summary>
/// Returns the lower Cholesky Factor, L, of input matrix A.
/// Satisfies the equation: L*L^T = A.
/// </summary>
/// <param name="a">Input matrix must be square, symmetric,
/// and positive definite. This method does not check for these properties,
/// and may produce unexpected results of those properties are not met.</param>
/// <returns></returns>
public static double[,] Cholesky(double[,] a)
{
int n = (int)Math.Sqrt(a.Length);
double[,] ret = new double[n, n];
for (int r = 0; r < n; r++)
for (int c = 0; c <= r; c++)
{
if (c == r)
{
double sum = 0;
for (int j = 0; j < c; j++)
{
sum += ret[c, j] * ret[c, j];
}
ret[c, c] = Math.Sqrt(a[c, c] - sum);
}
else
{
double sum = 0;
for (int j = 0; j < c; j++)
sum += ret[r, j] * ret[c, j];
ret[r, c] = 1.0 / ret[c, c] * (a[r, c] - sum);
}
}
return ret;
}
}
}
You may also check:How to resolve the algorithm Mandelbrot set step by step in the M4 programming language
You may also check:How to resolve the algorithm Faulhaber's triangle step by step in the Raku programming language
You may also check:How to resolve the algorithm Dijkstra's algorithm step by step in the M2000 Interpreter programming language
You may also check:How to resolve the algorithm Write float arrays to a text file step by step in the Lua programming language
You may also check:How to resolve the algorithm Hailstone sequence step by step in the Nim programming language