Skip to content

Instantly share code, notes, and snippets.

@pmslavin

pmslavin/cholesky.cpp

Last active June 28, 2022 22:14
Show Gist options
  • Save pmslavin/98fc5c03dd908bb4c7792992cd00ffd1 to your computer and use it in GitHub Desktop.
Save pmslavin/98fc5c03dd908bb4c7792992cd00ffd1 to your computer and use it in GitHub Desktop.
Download ZIP
Quadratic regression by Cholesky
Raw
cholesky.cpp
#include <vector>
#include <iostream>
#include <iomanip>
#include <numeric>
#include <cmath>
static const int PTERMS = 3;
/* Construct b vector of Ax=b */
double *build_bproduct(const std::vector<double>& xs, const std::vector<double>& ys)
{
int idx, n = xs.size();
double sy1 = std::accumulate(ys.begin(), ys.end(), 0.0);
double sx1y1=0.0, sx2y1=0.0;
for(idx=0; idx < n; ++idx){
sx1y1 += xs[idx]*ys[idx];
sx2y1 += xs[idx]*xs[idx]*ys[idx];
}
double *b = new double[PTERMS];
b[0] = sy1;
b[1] = sx1y1;
b[2] = sx2y1;
return b;
}
/* Construct coefficient array for degree 2 linear system */
double **build_2deg_coefficient_array(const std::vector<double>& xs)
{
int idx;
double n = double(xs.size());
double sx1 = std::accumulate(xs.begin(), xs.end(), 0.0);
double sx2=0.0, sx3=0.0, sx4=0.0;
double xsq;
for(auto x: xs){
xsq = x*x;
sx2 += xsq;
sx3 += xsq*x;
sx4 += xsq*xsq;
}
double **A = new double*[PTERMS];
for(idx=0; idx < PTERMS; ++idx)
A[idx] = new double[PTERMS];
A[0][0] = n;
A[0][1] = sx1;
A[0][2] = sx2;
A[1][0] = sx1;
A[1][1] = sx2;
A[1][2] = sx3;
A[2][0] = sx2;
A[2][1] = sx3;
A[2][2] = sx4;
return A;
}
/* Factors A into LL^ */
double **cholesky(double **A)
{
double **L = new double*[PTERMS];
for(int idx = 0; idx < PTERMS; ++idx)
L[idx] = new double[PTERMS];
L[0][0] = std::sqrt(A[0][0]);
L[0][1] = 0.0;
L[0][2] = 0.0;
L[1][0] = A[0][1]/L[0][0];
L[1][1] = std::sqrt(A[1][1] - L[1][0]*L[1][0]);
L[1][2] = 0.0;
L[2][0] = A[0][2]/L[0][0];
L[2][1] = (A[2][1] - L[2][0]*L[1][0])/L[1][1];
L[2][2] = std::sqrt(A[2][2] - (L[2][0]*L[2][0] + L[2][1]*L[2][1]));
return L;
}
/* Solve Ly=b via forward substitution */
double *solve_Lyb(double **L, double *b, std::vector<double>& ys)
{
double *y = new double[PTERMS];
double sy1 = std::accumulate(ys.begin(), ys.end(), 0.0);
y[0] = sy1/L[0][0];
y[1] = (b[1] - L[1][0]*y[0])/L[1][1];
y[2] = (b[2] - L[2][0]*y[0] - L[2][1]*y[1])/L[2][2];
return y;
}
/* Solve LTx=y via backwards substitution */
double *solve_LTxy(double **LT, double *y)
{
double *x = new double[PTERMS];
x[2] = y[2]/LT[2][2];
x[1] = (y[1] - LT[1][2]*x[2])/LT[1][1];
x[0] = (y[0] - LT[0][2]*x[2] - LT[0][1]*x[1])/LT[0][0];
return x;
}
/* Transpose square matrix L */
double **transpose(double **L)
{
double **LT = new double*[PTERMS];
for(int idx=0; idx < PTERMS; ++idx)
LT[idx] = new double[PTERMS];
for(int c=0; c < PTERMS; ++c)
for(int r=0; r < PTERMS; ++r)
LT[c][r] = L[r][c];
return LT;
}
void print_array(double **A)
{
for(int c=0; c < PTERMS; ++c){
for(int r=0; r < PTERMS; ++r){
std::cout << std::setw(8) << A[c][r] << " ";
}
std::cout << std::endl;
}
std::cout << std::endl;
}
void print_vector(double *b)
{
for(int idx=0; idx < PTERMS; ++idx)
std::cout << std::setw(8) << b[idx] << std::endl;
std::cout << std::endl;
}
/* Construct fit via solve of Ax=b, returning x */
double *fit_2deg_poly(std::vector<double>& xs, std::vector<double>& ys)
{
double **A = build_2deg_coefficient_array(xs);
double *b = build_bproduct(xs, ys);
print_array(A);
print_vector(b);
double **L = cholesky(A);
#ifdef DEBUG
print_array(L);
#endif
for(int idx=0; idx<PTERMS; ++idx)
delete [] A[idx];
delete [] A;
double *y = solve_Lyb(L, b, ys);
#ifdef DEBUG
print_vector(y);
#endif
delete [] b;
double **LT = transpose(L);
#ifdef DEBUG
print_array(LT);
#endif
for(int idx=0; idx<PTERMS; ++idx)
delete [] L[idx];
delete [] L;
double *x = solve_LTxy(LT, y);
for(int idx=0; idx<PTERMS; ++idx)
delete [] LT[idx];
delete [] LT;
delete [] y;
return x;
}
int main(void)
{
/* Test series: y=x^2 , y'=2x*/
// std::vector<double> xseries = { 0.0, 0.5, 1.0, 1.5, 2.0, 2.5 };
// std::vector<double> yseries = { 0.0, 0.25, 1.0, 2.25, 4.0, 6.25};
std::vector<double> xseries = { 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5 };
std::vector<double> yseries = { 11.0, 4.0, 4.0, 6.0, 9.0, 12.0, 17.0, 18.0, 14.0, 18.0, 13.0, 18.0};
double *x = fit_2deg_poly(xseries, yseries);
print_vector(x);
/* Fit func in std form */
auto f = [&x](double a){return x[0] + x[1]*a + x[2]*a*a;};
/* First deriv of fit func */
auto f1 = [&x](double a){return x[1] + 2.0*x[2]*a;};
/* Extrapolate to future t vals */
std::vector<double> tplus = {6.0, 6.5, 7.0};
for(auto& t: tplus){
std::cout << "f(" << std::fixed << std::setprecision(2) << t << ") = " \
<< std::setprecision(5) << std::setw(8) << f(t);
std::cout << std::setw(12) << std::fixed << std::setprecision(2) << "f'(" << std::setw(3) << t \
<< ") = " << std::setw(8) << std::setprecision(5) << f1(t) << std::endl;
}
delete [] x;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment