From 31eb41c51317e713b17ecf2b4fa2e05a35f00608 Mon Sep 17 00:00:00 2001
From: novak_99 <cooljohn312@gmail.com>
Date: Thu, 11 Nov 2021 21:48:49 -0800
Subject: [PATCH] cholesky decompositon! + positive semi-def checker for
 matrices

---
 MLPP/LinAlg/LinAlg.cpp | 39 ++++++++++++++++++++++++++++++++++++++-
 MLPP/LinAlg/LinAlg.hpp |  4 ++++
 main.cpp               | 28 +++++++++++++++++++++-------
 3 files changed, 63 insertions(+), 8 deletions(-)

diff --git a/MLPP/LinAlg/LinAlg.cpp b/MLPP/LinAlg/LinAlg.cpp
index 651cf71..0f22d8e 100644
--- a/MLPP/LinAlg/LinAlg.cpp
+++ b/MLPP/LinAlg/LinAlg.cpp
@@ -567,7 +567,7 @@ namespace MLPP{
 
         std::vector<std::vector<double>> a_new_prior = a_new;
         
-        // Bubble Sort
+        // Bubble Sort. Should change this later.
         for(int i = 0; i < a_new.size() - 1; i++){
             for(int j = 0; j < a_new.size() - 1 - i; j++){
                 if(a_new[j][j] < a_new[j + 1][j + 1]){
@@ -642,6 +642,29 @@ namespace MLPP{
         return {Q, R};
 
     }
+    
+    std::tuple<std::vector<std::vector<double>>, std::vector<std::vector<double>>> LinAlg::chol(std::vector<std::vector<double>> A){
+        std::vector<std::vector<double>> L = zeromat(A.size(), A[0].size());
+        for(int j = 0; j < L.size(); j++){ // Matrices entered must be square. No problem here.
+            for(int i = j; i < L.size(); i++){
+                if(i == j){
+                    double sum = 0;
+                    for(int k = 0; k < j; k++){ 
+                        sum += L[i][k] * L[i][k];
+                    }
+                    L[i][j] = std::sqrt(A[i][j] - sum);
+                }
+                else{ // That is, i!=j
+                    double sum = 0;
+                    for(int k = 0; k < j; k++){ 
+                        sum += L[i][k] * L[j][k];
+                    }
+                    L[i][j] = (A[i][j] - sum)/L[j][j];
+                }
+            }
+        }
+        return {L, transpose(L)}; // Indeed, L.T is our upper triangular matrix. 
+    }
 
     double LinAlg::sum_elements(std::vector<std::vector<double>> A){
         double sum = 0;
@@ -667,6 +690,20 @@ namespace MLPP{
         return mat_vec_mult(inverse(A), b);
     }
 
+    bool LinAlg::positiveDefiniteChecker(std::vector<std::vector<double>> A){
+        auto [eigenvectors, eigenvals] = eig(A);
+        std::vector<double> eigenvals_vec;
+        for(int i = 0; i < eigenvals.size(); i++){
+            eigenvals_vec.push_back(eigenvals[i][i]);
+        }
+        for(int i = 0; i < eigenvals_vec.size(); i++){
+            if(eigenvals_vec[i] <= 0){ // Simply check to ensure all eigenvalues are positive.
+                return false;
+            }
+        }
+        return true;
+    }
+
     void LinAlg::printMatrix(std::vector<std::vector<double>> A){
         for(int i = 0; i < A.size(); i++){
             for(int j = 0; j < A[i].size(); j++){
diff --git a/MLPP/LinAlg/LinAlg.hpp b/MLPP/LinAlg/LinAlg.hpp
index 410fd8a..c7a2992 100644
--- a/MLPP/LinAlg/LinAlg.hpp
+++ b/MLPP/LinAlg/LinAlg.hpp
@@ -94,11 +94,15 @@ namespace MLPP{
 
         std::tuple<std::vector<std::vector<double>>, std::vector<std::vector<double>>> QRD(std::vector<std::vector<double>> A);
 
+        std::tuple<std::vector<std::vector<double>>, std::vector<std::vector<double>>> chol(std::vector<std::vector<double>> A);
+
         double sum_elements(std::vector<std::vector<double>> A);
 
         std::vector<double> flatten(std::vector<std::vector<double>> A);
 
         std::vector<double> solve(std::vector<std::vector<double>> A, std::vector<double> b);
+
+        bool positiveDefiniteChecker(std::vector<std::vector<double>> A);
         
         void printMatrix(std::vector<std::vector<double>> A);
         
diff --git a/main.cpp b/main.cpp
index c36ad2b..870fc0d 100644
--- a/main.cpp
+++ b/main.cpp
@@ -433,18 +433,32 @@ int main() {
     // data.getImage("../../Data/apple.jpeg", chicken);
     // alg.printVector(chicken);
 
-    // TESTING QR DECOMP. EXAMPLE VIA WIKIPEDIA. SEE https://en.wikipedia.org/wiki/QR_decomposition.
+    // // TESTING QR DECOMP. EXAMPLE VIA WIKIPEDIA. SEE https://en.wikipedia.org/wiki/QR_decomposition.
 
-    std::vector<std::vector<double>> P = {{12, -51, 4}, {6, 167, -68}, {-4, 24, -41}};
-    alg.printMatrix(P);
+    // std::vector<std::vector<double>> P = {{12, -51, 4}, {6, 167, -68}, {-4, 24, -41}};
+    // alg.printMatrix(P);
 
-    alg.printMatrix(alg.gramSchmidtProcess(P));
+    // alg.printMatrix(alg.gramSchmidtProcess(P));
 
-    auto [Q, R] = alg.QRD(P); // It works! 
+    // auto [Q, R] = alg.QRD(P); // It works! 
     
-     alg.printMatrix(Q);
+    //  alg.printMatrix(Q);
 
-     alg.printMatrix(R); 
+    //  alg.printMatrix(R); 
+
+    // // Checking positive-definiteness checker. For Cholesky Decomp. 
+    // std::vector<std::vector<double>> A = 
+    // {
+    //     {1,-1,-1,-1},                        
+    //     {-1,2,2,2},
+    //     {-1,2,3,1},
+    //     {-1,2,1,4}
+    // };
+
+    // std::cout << std::boolalpha << alg.positiveDefiniteChecker(A) << std::endl;
+    // auto [L, Lt] = alg.chol(A); // WORKS !!!!
+    // alg.printMatrix(L);
+    // alg.printMatrix(Lt);
 
     return 0;
 }
\ No newline at end of file