用矩陣運算實現最小二乘法曲線擬合算法

1. 多項式擬合函數：

y= a0 + a1x + a2x^2 + ... + akx^k    (其中k為擬合次數)

當k=1 為線性擬合 ，k=2 為二次多項式 ...  三次多項式。

2. 最小二乘原理矩陣算法原理：

X*A=Y
A=((X'*X)-1)*X'*Y

|1  X1  X1^2  ... X1^k|            |y0|
|1  X2  X2^2  ... X2^k| |a0|     |y1|
|...                              | |a1|     |.  |
|...                              | |.   |  = |.  |
|...                              | |ak|     |.  |
|1  Xn  Xn^2  ... Xn^k|            |yn|

其中 X為初始范德蒙矩陣，A 為系數向量， Y 為因變量值向量

3.計算單元算法：幾個矩陣運算算法函數

1 ) 矩陣的轉置運算； 2)矩陣的逆運算；3 )矩陣乘法運算

4. C++代碼實現

1)矩陣乘法運算函數：

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）矩陣轉置運算函數：
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)矩陣逆運算函數：

// 工具函數

long double CMatrix::Det(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;
}
 

//工具函數
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];
        }
    }
 
}
 
// 逆運算函數
BOOL CMatrix::Inv(Matrix &a,Matrix &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）多項式擬合函數可以根據以上運算單元，結合矩陣運算公式：A=((X'*X)-1)*X'*Y

自由實現。

5）數據結構定義：

struct Matrix{  
    int m, n;   
    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;
        }
    }
};