用矩阵运算实现最小二乘法曲线拟合算法

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