最小二乘法

// test_2.cpp : Defines the entry point for the console application.
//最小二乘法

#include "stdafx.h"
#include
#include

const int n = 2;
const int m = 5;
int sgn(double x);
void LSS(double G[][n+1], int xm, int xn, double x[m], double &p, double w[m]);//LSS函数内部会修改G,x和p的值
                                                                                //该函数计算系数矩阵a[n]及误差p，系数存储在x[m]前n项中

int _tmain(int argc, _TCHAR* argv[])
{
    using namespace std;
    double temp = 0.0;
    double p = 0.0;//保存误差
    double G[m][n+1] = {0};
    double x[m] = {0.1, 0.2, 0.3, 0.4, 0.5};//x[n]同时存储最终结果a[n]
    double w[m] = {0.0};

    //给G[m][n+1]赋值
    for(int i = 0; i < m; i++)
    {
        G[i][0] = 1.0;
        G[i][1] = x[i];
    }
    //将y值写入G[m][n+1]最后一列
    G[0][n] = 5.1234;
    G[1][n] = 5.3057;
    G[2][n] = 5.5678;
    G[3][n] = 5.9378;
    G[4][n] = 6.4370;
    //调用最小二乘函数
    LSS(G, m, n, x, p, w);

    //4 输出结果
    //输出系数
    cout << "系数为：/n";
    for(int i = 0; i < n; i++)
    {
        cout << "a" << i << " = " << x[i] << endl;
    }

    //输出误差
    cout << "/n误差为 " << p << endl;

    return 0;
}

 

void LSS(double G[][n+1], int xm, int xn, double x[m], double &p, double w[m])
{
    double b = 0.0;//保存二范数
    double B = 0.0;
    double t = 0.0;
    double temp = 0.0;
    //double w[m] = {0.0};//这里不能是xm，这样不好！！！

    //最小二乘算法
    //1
    for(int k = 0; k < xn; k++)
    {
        //1.1 形成矩阵Qk
        temp = 0.0;
        for(int i = k; i < xm; i++)
            temp += G[i][k] * G[i][k];
        b = -sgn(G[k][k]) * sqrt(temp);
        w[k] = G[k][k] - b;
        for(int j = k+1; j < xm; j++)
            w[j] = G[j][k];
        B = b * w[k];
        //1.2 变换Gk-1到Gk
        G[k][k] = b;
        for(int j = k+1; j < xn+1; j++)
        {   
            temp = 0.0;
            for(int i = k; i < xm; i++)
                temp += w[i] * G[i][j];
            t = temp / B;
            for(int i = k; i < xm; i++)
            G[i][j] += t * w[i];
        }
    }

    //2 解三角方程Rx=h1
    x[xn-1] = G[xn-1][xn] / G[xn-1][xn-1];
    for(int i = xn-1; i >= 0; i--)
    {
        temp = 0.0;
        for(int j = i+1; j < xn; j++)
            temp += G[i][j] * x[j];
        x[i] = (G[i][xn] - temp) / G[i][i];
    }

    //3 计算误差E2
    for(int i = xn; i < xm; i++)
        p += G[i][xn] * G[i][xn];
}

 

int sgn(double x)
{
    if (x>0) return(1);    //返回出口1:如果键盘输入正数，则返回1
    if (x<0) return(-1);   //返回出口2:如果键盘输入负数，则返回-1
    return(0);             //返回出口3:如果键盘输入0，   则返回0
}

 

最小二乘法的算法原理见《计算方法教程（第二版）》 凌永祥、陈明奎 西安交通大学出版社 2005年4月第2版 P123