Introduction:
Cholesky Factorization is a computational method commonly applied in linear algebra to break down a Hermitian, positive-definite matrix into the multiplication of a lower triangular matrix and its conjugate transpose. This approach is notably effective in resolving linear equation systems, calculating determinants, and executing numerical simulations.
In C++, executing Cholesky Decomposition entails decomposing the matrix into lower triangular parts, streamlining future computations and enhancements. The process involves iteratively computing the Cholesky factor by adjusting the elements within the lower triangular matrix.
Implementing Cholesky Decomposition in C++ often requires creating classes or functions that can execute matrix operations effectively, such as matrix multiplication, transposition, and solving linear systems. It is essential to also pay attention to error handling and memory allocation for a well-rounded implementation.
Program:
#include <iostream>
#include <vector>
#include <cmath>
// Function to perform Cholesky Decomposition
std::vector<std::vector<double>> choleskyDecomposition(const std::vector<std::vector<double>>& matrix) {
int n = matrix.size();
std::vector<std::vector<double>> L(n, std::vector<double>(n, 0.0));
for (int i = 0; i < n; ++i) {
for (int j = 0; j <= i; ++j) {
double sum = 0.0;
if (j == i) {
for (int k = 0; k < j; ++k) {
sum += pow(L[j][k], 2);
}
L[j][j] = sqrt(matrix[j][j] - sum);
} else {
for (int k = 0; k < j; ++k) {
sum += (L[i][k] * L[j][k]);
}
L[i][j] = (matrix[i][j] - sum) / L[j][j];
}
}
}
return L;
}
// Function to display a matrix
void displayMatrix(const std::vector<std::vector<double>>& matrix) {
for (const auto& row : matrix) {
for (double element : row) {
std::cout << element << " ";
}
std::cout << std::endl;
}
}
int main() {
// Example matrix
std::vector<std::vector<double>> A = {{4, 12, -16},
{12, 37, -43},
{-16, -43, 98}};
std::cout << "Original Matrix A:" << std::endl;
displayMatrix(A);
// Perform Cholesky Decomposition
std::vector<std::vector<double>> L = choleskyDecomposition(A);
std::cout << "\nCholesky Decomposition L:" << std::endl;
displayMatrix(L);
return 0;
}
Output:
Original Matrix A:
4 12 -16
12 37 -43
-16 -43 98
Cholesky Decomposition L:
2 0 0
6 1 0
-8 5 3
Explanation:
- Matrix Representation:
The matrix is depicted as a two-dimensional array of double values (std::vector<std::vector<double>>). This setup enables convenient handling of the matrix components.
- Cholesky Decomposition Method (choleskyDecomposition):
The function accepts a matrix as an argument and produces the lower triangular Cholesky factor matrix as output.
It creates a new matrix L initialized with zeros to save the Cholesky factor.
It loops through every row and column within the matrix, computing the values of the Cholesky decomposition matrix L.
At every step, it calculates the values in L using the Cholesky factorization method.
The algorithm computes every entry in the lower triangular matrix L by adding the elements from the original matrix and the previously derived entries of L.
If the existing element falls on the main diagonal, its value is determined by taking the square root of the disparity between the corresponding element in the initial matrix and the total of squares of previously computed L elements.
If the current item is not located along the diagonal, it is computed by utilizing the values of L that were calculated earlier.
- Implement the Display Matrix Function (displayMatrix):
This function is utilized to showcase the contents of a matrix on the console.
It loops over every row in the matrix and displays each item with spaces between them.
- Main Function:
It defines an example matrix A.
It displays the original matrix A.
It invokes the choleskyDecomposition function to execute Cholesky Decomposition on matrix A.
It displays the Cholesky factor matrix L.
Complexity analysis:
Time Complexity:
- The main operations in the Cholesky Decomposition algorithm involve nested loops iterating over the rows and columns of the matrix.
- The outer loop runs n times, where n is the number of rows or columns in the matrix.
- The inner loops run at most i times for each iteration of the outer loop, where i is the index of the current row/column.
- Within the inner loops, there are computations involving matrix multiplications and additions, which are typically constant time operations.
- Therefore, the overall time complexity of the Cholesky Decomposition algorithm is approximately O(n3) , where n is the size of the matrix.
Space Complexity:
- The space complexity of the Cholesky Decomposition algorithm mainly depends on the space required to store the Cholesky factor matrix L.
- The Cholesky factor matrix L is a lower triangular matrix of size n×n.
- Since only the lower triangular part of the matrix is stored, the space complexity is approximately O(n2) .
- Additionally, there is some auxiliary space required for variables such as loop counters and temporary storage, but this is negligible compared to the space required for the Cholesky factor matrix.
- Therefore, the overall space complexity of the Cholesky Decomposition algorithm is approximately O(n2) .
Blocked Cholesky Decomposition
Blocked Cholesky Decomposition divides the matrix into smaller blocks and applies the Cholesky Decomposition to each block, leveraging cache proximity and minimizing memory access operations.
This method consists of partitioning the initial matrix into square submatrices, also known as blocks, and conducting Cholesky Decomposition on each block separately.
By breaking down smaller blocks, the algorithm can effectively utilize processor caches and enhance memory access efficiency. Once the blocks are decomposed, the lower triangular sections are combined to create the Cholesky factor matrix.
Program:
#include <iostream>
#include <vector>
#include <cmath>
// Function to perform Cholesky Decomposition on a block
void choleskyBlock(std::vector<std::vector<double>>& A, int start, int end) {
for (int j = start; j < end; ++j) {
for (int k = 0; k < j; ++k) {
for (int i = j; i < end; ++i) {
A[i][j] -= A[i][k] * A[j][k];
}
}
A[j][j] = sqrt(A[j][j]);
for (int i = j + 1; i < end; ++i) {
A[i][j] /= A[j][j];
}
}
}
// Function to perform Blocked Cholesky Decomposition
void choleskyBlocked(std::vector<std::vector<double>>& A, int blockSize) {
int n = A.size();
for (int k = 0; k < n; k += blockSize) {
int end = std::min(k + blockSize, n);
choleskyBlock(A, k, end);
for (int i = end; i < n; ++i) {
for (int j = k; j < end; ++j) {
for (int p = 0; p < end; ++p) {
A[i][j] -= A[i][p] * A[j][p];
}
}
}
}
}
// Function to display a matrix
void displayMatrix(const std::vector<std::vector<double>>& matrix) {
for (const auto& row : matrix) {
for (double element : row) {
std::cout << element << " ";
}
std::cout << std::endl;
}
}
int main() {
// Example matrix
std::vector<std::vector<double>> A = {{25, 15, -5},
{15, 18, 0},
{-5, 0, 11}};
int blockSize = 2; // Block size for blocked Cholesky Decomposition
std::cout << "Original Matrix A:" << std::endl;
displayMatrix(A);
// Perform Blocked Cholesky Decomposition
choleskyBlocked(A, blockSize);
std::cout << "\nCholesky Decomposition L:" << std::endl;
displayMatrix(A);
return 0;
}
Output:
Original Matrix A:
25 15 -5
15 18 0
-5 0 11
Cholesky Decomposition L:
5 15 -5
3 3 0
20 120 -nan
Explanation:
- choleskyBlock Function:
It executes Cholesky Decomposition on a section of the matrix.
It loops over the rows and columns within the block, modifying each element according to the Cholesky Decomposition method.
Elements undergo modifications within nested loops, where the operations of multiplication and subtraction are managed.
- choleskyBlocked Function:
Splits the matrix into smaller sections and conducts Cholesky Decomposition on each individual section.
Traverses the matrix in block-sized increments and invokes the choleskyBlock function for each block.
Manages the final block, which might be smaller in size compared to the designated block size.
- displayMatrix Function:
Displays the contents of a matrix to the console.
Goes through every row in the matrix and displays each item with a space in between.
- Main Function:
Defines an example matrix A.
Specifies the size of blocks for performing Blocked Cholesky Decomposition.
Displays the original matrix A.
Invokes the choleskyBlocked function to execute Blocked Cholesky Decomposition on matrix A.
Displays the resulting Cholesky factor matrix.
Complexity analysis:
Time Complexity:
- The main operations in the algorithm involve iterating over the matrix elements and performing arithmetic operations within nested loops.
- The outer loop in the choleskyBlocked function runs O(n/blockSize) times, where n is the size of the matrix and blockSize is the specified block size.
- Within each block, the choleskyBlock function performs Cholesky Decomposition, which involves additional nested loops. The number of iterations for each block depends on the block size and is relatively constant.
- Overall, the time complexity of the Blocked Cholesky Decomposition algorithm is approximately O(n2/blockSize) .
Space Complexity:
- The space complexity primarily depends on the storage of the input matrix and the Cholesky factor matrix.
- Since the input matrix A and the Cholesky factor matrix are both represented as two-dimensional vectors of doubles, their combined space complexity is O(n2) , where n is the size of the matrix.
- Additionally, there is some auxiliary space required for loop counters and temporary variables within functions, but this is negligible compared to the size of the matrices.
- Therefore, the overall space complexity of the Blocked Cholesky Decomposition algorithm is approximately O(n2) .
