Using matrix operation to realize least squares curve fitting algorithm

1. Polynomial-fitting function:

Y= A0 + a1x + a2x^2 + ... + akx^k (where k is fit times)

When the k=1 is linear fit, the k=2 is a two-time polynomial ... Three-time polynomial.

2. Least squares principle matrix algorithm principle:

X*a=y
A= ((X ' *x)-1) *x ' *y

|-X1 x1^2 ...            x1^k| |y0|
|-X2 x2^2 ...     x2^k| |a0| |y1|
|...                              |     |a1|  |. |
|...                              |   |.  |  = |. |
|...                              |     |ak|  |. |
|-Xn xn^2 ...            xn^k| |yn|

where x is the initial Vandermonde matrix, A is the coefficient vector, Y is the dependent variable value vector

3. Computational unit algorithm: several matrix arithmetic functions

1) Matrix's transpose operation; 2) matrix inverse; 3) matrix multiplication

4. C + + code implementation

1) matrix multiplication operation function:

BOOL Cmatrix::mul (const matrix &a, const matrix &b, Matrix &c)
{
if (A.N!=B.M)
{
return FALSE;
}

C.M = a.m.;
C.N = B.N;
int I, j, K;
for (i=0; i<a.m; i++)
{
for (j=0; j<b.n; j + +)
{
For (C.p[i * C.N + j] = 0, K =0; k<a.n; k++)
{
C.p[i * C.N + j] + = A.p[i * A.N + K] * B.p[k * B.N + j];
}
}
}
return TRUE;
}
2) Matrix transpose operation function:
void Cmatrix::trs (Matrix &a, Matrix &trs)
{
TRS.M = A.N;
TRS.N = a.m.;

for (int i = 0; i < a.m.; i++)
{
for (int j = 0; J < A.N; J + +)
{
Trs.p[j * a.m. + i] = A.p[i * A.N + j];
}
}
}

3) matrix inverse function:

Tool functions

Long double Cmatrix::D et (Matrix &a)
{
Long double det = 0;
if (a.m.! = A.N)
{
//...
return 0;
}
if (A.N = = 1)
{
Det = a.p[0];
return det;
}
Else
{
for (int i = 0; i < A.N; i++)
{
if (i% 2 = = 0)
{
Matrix Mat;
Adjunct (A, I, 0,mat);
Det + = a.p[i * A.N] * DET (MAT);
Delete MAT.P;
}
Else
{
Matrix Mat;
Adjunct (A, I, 0,mat);
Det-= a.p[i * A.N] * DET (MAT);
Delete MAT.P;
}
}
}

return det;
}

Tool functions
void Cmatrix::adjunct (Matrix A, int indexm, int indexn,matrix &adj)
{
Adj. SetSize (a.n-1, a.n-1);
for (int i = 0; i < INDEXM; i++)
{
for (int j = 0; J < Indexn; J + +)
{
Adj.p[i * (a.n-1) + j] = A.p[i * A.N + j];
}
for (int k = indexn + 1; k < A.N; k++)
{
Adj.p[i * (a.n-1) + k-1] = A.p[i * A.N + K];
}
}

for (int m = INDEXM + 1; m < A.N; m++)
{
for (int j = 0; J < A.n-1; J + +)
{
adj.p[(m-1) * (a.n-1) + j] = A.p[m * A.N + j];
}
for (int k = indexn + 1; k < A.N; k++)
{
adj.p[(m-1) * (a.n-1) + k-1] = A.p[m * A.N + K];
}
}

}

Inverse function
BOOL CMATRIX::INV (Matrix &a,matrix &AMP;INV)
{
Matrix Temp (A.M,A.N);

if (A.M!=A.N)
{
return FALSE;
}


Long double det = det (a);
if (Det = = 0)
{
return FALSE;
}

for (int i = 0; i < temp.m; i++)
{
for (int j = 0; J < TEMP.N; J + +)
{
if ((i +j)% 2 = = 0)
{
Matrix Mat;
Adjunct (A, I, J,mat);
Temp.p[i * temp.m + j] = Det (MAT)/Det;
Delete MAT.P;
}
Else
{
Matrix Mat;
Adjunct (A, I, J,mat);
Temp.p[i * temp.m + j] =-det (MAT)/Det;
Delete MAT.P;
}
}
}

Trs (TEMP,INV);
return TRUE;
}

4) Polynomial fitting function can be combined with matrix operation formula according to the above arithmetic unit:a= ((X ' *x)-1) *x ' *y

Free to achieve.

5) Data structure definition:

struct matrix{
int m, n;
a long double *p;
Matrix ()
{
m = 0;
n = 0;
p = NULL;
}
Matrix (int t_m, int t_n)
{
m = t_m;
n = t_n;
p = new long double[m*n];
}
void SetSize (int t_m,int t_n)
{
m = t_m;
n = t_n;
p = new long double[m*n];
}
~matrix ()
{
if (P! = NULL)
{
Delete p;
}
}
};

Using matrix operation to realize least squares curve fitting algorithm