cholesky decompositon! + positive semi-def checker for matrices

This commit is contained in:
novak_99 2021-11-11 21:48:49 -08:00
parent 8687fc5eb1
commit 31eb41c513
3 changed files with 63 additions and 8 deletions

View File

@ -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++){

View File

@ -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);

View File

@ -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;
}