最小二乘法用于直线，多项式，圆，椭圆的拟合及程序实现

最小二乘法是一种优化算法，最小二乘法名字的缘由有两个：一是要将误差最小化，二是将误差最小化的方法是使误差的平方和最小化。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合，所拟合的曲线可以是线性拟合与非线性拟合。

--------------------一元线性函数--------------------

形式1：借用案例（http://blog.csdn.net/qll125596718/article/details/8248249）首先以一元线性方程参数估计为例，样本回归模型：

残差平方和：

则通过Q最小确定这条直线，即确定，以为变量，把它们看作是Q的函数，就变成了一个求极值的问题，可以通过求导数得到。求Q对两个待估参数的偏导数：

解得：

实例（c++）：

#include <iostream>
#include <algorithm>
#include <valarray>
#include <vector>
using namespace std;int main()
{double x[] = { 1, 2, 3, 4, 5, 6 };double y[] = { 3, 5.5, 6.8, 8.8, 11, 12};valarray<double> data_x(x, 6);valarray<double> data_y(y, 6);float A = 0.0;float B = 0.0;float C = 0.0;float D = 0.0;A = (data_x*data_x).sum();B = data_x.sum();C = (data_x*data_y).sum();D = data_y.sum();float tmp = A*data_x.size() - B*B;float k, b;k = (C*data_x.size() - B*D) / tmp;b = (A*D - C*B) / tmp;cout << "y=" << k << "x+" << b << endl;return 0;
}

运行结果：

注：valarray类似vector，也是一个模板类，主要用来对一系列元素进行高速的数字计算，其与vector的主要区别在于以下两点：
1、valarray定义了一组在两个相同长度和相同类型的valarray类对象之间的数字计算；
2、通过重载operater[]，可以返回valarray的相关信息（valarray其中某个元素的引用、特定下标的值或者其某个子集）。

valarray类构造函数:
valarray( );
explicit valarray(size_t _Count);
valarray( const Type& _Val, size_t _Count);
valarray( const Type* _Ptr, size_t _Count);
valarray( const slice_array<Type>& _SliceArray);
valarray( const gslice_array<Type>& _GsliceArray);
valarray( const mask_array<Type>& _MaskArray);
valarray( const indirect_array<Type>& _IndArray);

valarray 类用法:
1. apply 将 valarray 数组的每一个值都用 apply 所接受到的函数进行计算
2. cshift 将 valarray 数组的数据进行循环移动，参数为正者左移为负就右移
3. max 返回 valarray 数组的最大值
4. min 返回 valarray 数组的最小值
5. resize 重新设置 valarray 数组大小，并对其进行初始化
6. shift 将 valarray 数组移动，参数为正者左移，为负者右移，移动后由 0 填充剩余位
7. size 得到数组的大小
8. sum 数组求和

--------------------N元线性函数--------------------

一元线性方程可以看做多元函数的特例，现在用矩阵形式表述多元函数情况下，最小二乘的一般形式：

设拟合多项式为：

各点到这条曲线的距离之和，即偏差平方和如下：

对等式右边求ai偏导数，得到：

.......

把这些等式表示成矩阵的形式，就可以得到下面的矩阵：

（3）

进行化简计算：

，

上面公式（3）可以写为:

，

#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "vector"
using namespace std;struct point
{double x;double y;
};typedef vector<double> doubleVector;
vector<point> getFile(char *File);  //获取文件数据
doubleVector getCoeff(vector<point> sample, int n);   //矩阵方程void main()
{int i, n;char *File = "XY.txt";vector<point> sample;doubleVector  coefficient;sample = getFile(File);printf("拟合多项式阶数n=");scanf_s("%d", &n);coefficient = getCoeff(sample, n);printf("\n拟合矩阵的系数为：\n");for (i = 0; i < coefficient.size(); i++)printf("a%d = %lf\n", i, coefficient[i]);
}
//矩阵方程
doubleVector getCoeff(vector<point> sample, int n)
{vector<doubleVector> matFunX;  //公式3左矩阵vector<doubleVector> matFunY;  //公式3右矩阵doubleVector temp;double sum;int i, j, k;//公式3左矩阵for (i = 0; i <= n; i++){temp.clear();for (j = 0; j <= n; j++){sum = 0;for (k = 0; k < sample.size(); k++)sum += pow(sample[k].x, j + i);temp.push_back(sum);}matFunX.push_back(temp);}//printf("matFunX.size=%d\n", matFunX.size());//printf("matFunX[3][3]=%f\n", matFunX[3][3]);//公式3右矩阵for (i = 0; i <= n; i++){temp.clear();sum = 0;for (k = 0; k < sample.size(); k++)sum += sample[k].y*pow(sample[k].x, i);temp.push_back(sum);matFunY.push_back(temp);}printf("matFunY.size=%d\n", matFunY.size());//矩阵行列式变换double num1, num2, ratio;for (i = 0; i < matFunX.size() - 1; i++){num1 = matFunX[i][i];for (j = i + 1; j < matFunX.size(); j++){num2 = matFunX[j][i];ratio = num2 / num1;for (k = 0; k < matFunX.size(); k++)matFunX[j][k] = matFunX[j][k] - matFunX[i][k] * ratio;matFunY[j][0] = matFunY[j][0] - matFunY[i][0] * ratio;}}//计算拟合曲线的系数doubleVector coeff(matFunX.size(), 0);for (i = matFunX.size() - 1; i >= 0; i--){if (i == matFunX.size() - 1)coeff[i] = matFunY[i][0] / matFunX[i][i];else{for (j = i + 1; j < matFunX.size(); j++)matFunY[i][0] = matFunY[i][0] - coeff[j] * matFunX[i][j];coeff[i] = matFunY[i][0] / matFunX[i][i];}}return coeff;
}//获取文件数据
vector<point> getFile(char *File)
{int i = 1;vector<point> dst;FILE *fp=fopen(File, "r");if (fp == NULL){printf("Open file error!!!\n");exit(0);}point temp;double num;while (fscanf(fp, "%lf", &num) != EOF){if (i % 2 == 0){temp.y = num;dst.push_back(temp);}elsetemp.x = num;i++;}fclose(fp);return dst;
}

XY.txt内容：

另外在http://blog.csdn.net/lsh_2013/article/details/46697625里也有相关程序：

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;  //最小二乘拟合相关函数定义
double sum(vector<double> Vnum, int n);
double MutilSum(vector<double> Vx, vector<double> Vy, int n);
double RelatePow(vector<double> Vx, int n, int ex);
double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex);
void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[]);
void CalEquation(int exp, double coefficient[]);
double F(double c[],int l,int m);
double Em[6][4];  //主函数，这里将数据拟合成二次曲线
int main(int argc, char* argv[])
{  double arry1[5]={0,0.25,0,5,0.75};  double arry2[5]={1,1.283,1.649,2.212,2.178};  double coefficient[5];  memset(coefficient,0,sizeof(double)*5);  vector<double> vx,vy;  for (int i=0; i<5; i++)  {  vx.push_back(arry1[i]);  vy.push_back(arry2[i]);  }  EMatrix(vx,vy,5,3,coefficient);  printf("拟合方程为：y = %lf + %lfx + %lfx^2 \n",coefficient[1],coefficient[2],coefficient[3]);  return 0;
}
//累加
double sum(vector<double> Vnum, int n)
{  double dsum=0;  for (int i=0; i<n; i++)  {  dsum+=Vnum[i];  }  return dsum;
}
//乘积和
double MutilSum(vector<double> Vx, vector<double> Vy, int n)
{  double dMultiSum=0;  for (int i=0; i<n; i++)  {  dMultiSum+=Vx[i]*Vy[i];  }  return dMultiSum;
}
//ex次方和
double RelatePow(vector<double> Vx, int n, int ex)
{  double ReSum=0;  for (int i=0; i<n; i++)  {  ReSum+=pow(Vx[i],ex);  }  return ReSum;
}
//x的ex次方与y的乘积的累加
double RelateMutiXY(vector<double> Vx, vector<double> Vy, int n, int ex)
{  double dReMultiSum=0;  for (int i=0; i<n; i++)  {  dReMultiSum+=pow(Vx[i],ex)*Vy[i];  }  return dReMultiSum;
}
//计算方程组的增广矩阵
void EMatrix(vector<double> Vx, vector<double> Vy, int n, int ex, double coefficient[])
{  for (int i=1; i<=ex; i++)  {  for (int j=1; j<=ex; j++)  {  Em[i][j]=RelatePow(Vx,n,i+j-2);  }  Em[i][ex+1]=RelateMutiXY(Vx,Vy,n,i-1);  }  Em[1][1]=n;  CalEquation(ex,coefficient);
}
//求解方程
void CalEquation(int exp, double coefficient[])
{  for(int k=1;k<exp;k++) //消元过程  {  for(int i=k+1;i<exp+1;i++)  {  double p1=0;  if(Em[k][k]!=0)  p1=Em[i][k]/Em[k][k];  for(int j=k;j<exp+2;j++)   Em[i][j]=Em[i][j]-Em[k][j]*p1;  }  }  coefficient[exp]=Em[exp][exp+1]/Em[exp][exp];  for(int l=exp-1;l>=1;l--)   //回代求解  coefficient[l]=(Em[l][exp+1]-F(coefficient,l+1,exp))/Em[l][l];
}
//供CalEquation函数调用
double F(double c[],int l,int m)
{  double sum=0;  for(int i=l;i<=m;i++)  sum+=Em[l-1][i]*c[i];  return sum;
}

memset相关介绍：

http://baike.baidu.com/link?url=p7JreiRCj9yPs3r3WAfsXgynjvtGrWoQ_exF9tFGK6fsVP7V6tdm-_13QhCZxqPrfRi0wH0EihhRL_-qVvrewq

http://c.biancheng.net/cpp/html/157.html

--------------------拟合圆的方程--------------------

/*
最小二乘法拟合圆，拟合出的圆以圆心坐标和半径的形式表示*/
typedef complex<int> POINT;
bool FitCircle(const std::vector<POINT> &points, double &er_x, double &er_y, double &radius)
{cent_x = 0.0f;cent_y = 0.0f;radius = 0.0f;if (points.size() < 3){return false;}double sum_x = 0.0f, sum_y = 0.0f;double sum_x2 = 0.0f, sum_y2 = 0.0f;double sum_x3 = 0.0f, sum_y3 = 0.0f;double sum_xy = 0.0f, sum_x1y2 = 0.0f, sum_x2y1 = 0.0f;int N = points.size();for (int i = 0; i < N; i++){double x = points[i].real();double y = points[i].imag();double x2 = x * x;double y2 = y * y;sum_x += x;sum_y += y;sum_x2 += x2;sum_y2 += y2;sum_x3 += x2 * x;sum_y3 += y2 * y;sum_xy += x * y;sum_x1y2 += x * y2;sum_x2y1 += x2 * y;}double C, D, E, G, H;double a, b, c;C = N * sum_x2 - sum_x * sum_x;D = N * sum_xy - sum_x * sum_y;E = N * sum_x3 + N * sum_x1y2 - (sum_x2 + sum_y2) * sum_x;G = N * sum_y2 - sum_y * sum_y;H = N * sum_x2y1 + N * sum_y3 - (sum_x2 + sum_y2) * sum_y;a = (H * D - E * G) / (C * G - D * D);b = (H * C - E * D) / (D * D - G * C);c = -(a * sum_x + b * sum_y + sum_x2 + sum_y2) / N;cent_x = a / (-2);cent_y = b / (-2);radius = sqrt(a * a + b * b - 4 * c) / 2;return true;
}

--------------------拟合椭圆方程--------------------

//LSEllipse.h

/*************************************************************************功能说明： 对平面上的一些列点给出最小二乘的椭圆拟合，利用奇异值分解法解得最小二乘解作为椭圆参数。调用形式： cvFitEllipse2f（arrayx,arrayy,box）；    参数说明： arrayx: arrayx[n],每个值为x轴一个点arrayx: arrayy[n],每个值为y轴一个点n     : 点的个数box   : box[5],椭圆的五个参数，分别为center.x,center.y,2a,2b,xthetaesp: 解精度，通常取1e-6,这个是解方程用的说
***************************************************************************/
#include<cstdlib>
#include<float.h>
#include<vector>
using namespace std;class LSEllipse
{
public:LSEllipse(void);~LSEllipse(void);vector<double> getEllipseparGauss(vector<CPoint> vec_point);void cvFitEllipse2f( int *arrayx, int *arrayy,int n,float *box );
private:int SVD(float *a,int m,int n,float b[],float x[],float esp);int gmiv(float a[],int m,int n,float b[],float x[],float aa[],float eps,float u[],float v[],int ka);int ginv(float a[],int m,int n,float aa[],float eps,float u[],float v[],int ka);int muav(float a[],int m,int n,float u[],float v[],float eps,int ka);
};

//LSEllipse.cpp

#include "LSEllipse.h"
#include <cmath>LSEllipse::LSEllipse(void)
{
}LSEllipse::~LSEllipse(void)
{
}
//列主元高斯消去法
//A为系数矩阵，x为解向量，若成功，返回true，否则返回false，并将x清空。bool RGauss(const vector<vector<double> > & A, vector<double> & x)
{x.clear();int n = A.size();int m = A[0].size();x.resize(n);//复制系数矩阵，防止修改原矩阵vector<vector<double> > Atemp(n);for (int i = 0; i < n; i++){vector<double> temp(m);for (int j = 0; j < m; j++){temp[j] = A[i][j];}Atemp[i] = temp;temp.clear();}for (int k = 0; k < n; k++){//选主元double max = -1;int l = -1;for (int i = k; i < n; i++){if (abs(Atemp[i][k]) > max){max = abs(Atemp[i][k]);l = i;}}if (l != k){//交换系数矩阵的l行和k行for (int i = 0; i < m; i++){double temp = Atemp[l][i];Atemp[l][i] = Atemp[k][i];Atemp[k][i] = temp;}}//消元for (int i = k+1; i < n; i++){double l = Atemp[i][k]/Atemp[k][k];for (int j = k; j < m; j++){Atemp[i][j] = Atemp[i][j] - l*Atemp[k][j];}}}//回代x[n-1] = Atemp[n-1][m-1]/Atemp[n-1][m-2];for (int k = n-2; k >= 0; k--){double s = 0.0;for (int j = k+1; j < n; j++){s += Atemp[k][j]*x[j];}x[k] = (Atemp[k][m-1] - s)/Atemp[k][k];}return true;
}vector<double>  LSEllipse::getEllipseparGauss(vector<CPoint> vec_point)
{vector<double> vec_result;double x3y1 = 0,x1y3= 0,x2y2= 0,yyy4= 0, xxx3= 0,xxx2= 0,x2y1= 0,yyy3= 0,x1y2= 0 ,yyy2= 0,x1y1= 0,xxx1= 0,yyy1= 0;int N = vec_point.size();for (int m_i = 0;m_i < N ;++m_i ){double xi = vec_point[m_i].x ;double yi = vec_point[m_i].y;x3y1 += xi*xi*xi*yi ;x1y3 += xi*yi*yi*yi;x2y2 += xi*xi*yi*yi; ;yyy4 +=yi*yi*yi*yi;xxx3 += xi*xi*xi ;xxx2 += xi*xi ;x2y1 += xi*xi*yi;x1y2 += xi*yi*yi;yyy2 += yi*yi;x1y1 += xi*yi;xxx1 += xi;yyy1 += yi;yyy3 += yi*yi*yi;}double resul[5];resul[0] = -(x3y1);resul[1] = -(x2y2);resul[2] = -(xxx3);resul[3] = -(x2y1);resul[4] = -(xxx2);long double Bb[5],Cc[5],Dd[5],Ee[5],Aa[5];Bb[0] = x1y3, Cc[0] = x2y1, Dd[0] = x1y2, Ee[0] = x1y1, Aa[0] = x2y2;Bb[1] = yyy4, Cc[1] = x1y2, Dd[1] = yyy3, Ee[1] = yyy2, Aa[1] = x1y3;Bb[2] = x1y2, Cc[2] = xxx2, Dd[2] = x1y1, Ee[2] = xxx1, Aa[2] = x2y1;Bb[3] = yyy3, Cc[3]= x1y1, Dd[3] = yyy2, Ee[3] = yyy1, Aa[3] = x1y2;Bb[4]= yyy2, Cc[4]= xxx1, Dd[4] = yyy1, Ee[4] = N, Aa[4]= x1y1;vector<vector<double>>Ma(5);vector<double>Md(5);for(int i=0;i<5;i++){Ma[i].push_back(Aa[i]);Ma[i].push_back(Bb[i]);Ma[i].push_back(Cc[i]);Ma[i].push_back(Dd[i]);Ma[i].push_back(Ee[i]);Ma[i].push_back(resul[i]);}RGauss(Ma,Md);long double A=Md[0];long double B=Md[1];long double C=Md[2];long double D=Md[3];long double E=Md[4];double XC=(2*B*C-A*D)/(A*A-4*B);double YC=(2*D-A*C)/(A*A-4*B);long double fenzi=2*(A*C*D-B*C*C-D*D+4*E*B-A*A*E);long double fenmu=(A*A-4*B)*(B-sqrt(A*A+(1-B)*(1-B))+1);long double fenmu2=(A*A-4*B)*(B+sqrt(A*A+(1-B)*(1-B))+1);double XA=sqrt(fabs(fenzi/fenmu));double XB=sqrt(fabs(fenzi/fenmu2));double Xtheta=0.5*atan(A/(1-B))*180/3.1415926;if(B<1)Xtheta+=90;vec_result.push_back(XC);vec_result.push_back(YC);vec_result.push_back(XA);vec_result.push_back(XB);vec_result.push_back(Xtheta);return vec_result;
}void  LSEllipse::cvFitEllipse2f(  int *arrayx, int *arrayy,int n,float *box )
{   float cx=0,cy=0;double rp[5], t;float *A1=new float[n*5];float *A2=new float[2*2];float *A3=new float[n*3];float *B1=new float[n],*B2=new float[2],*B3=new float[n];const double min_eps = 1e-6;int i;for( i = 0; i < n; i++ ){cx += arrayx[i]*1.0;cy += arrayy[i]*1.0;}cx /= n;cy /= n;for( i = 0; i < n; i++ ){int step=i*5;float px,py;px = arrayx[i]*1.0;py = arrayy[i]*1.0;px -= cx;py -= cy;B1[i] = 10000.0;A1[step] = -px * px;A1[step + 1] = -py * py;A1[step + 2] = -px * py;A1[step + 3] = px;A1[step + 4] = py;}float *x1=new float[5];//解出Ax^2+By^2+Cxy+Dx+Ey=10000的最小二乘解！SVD(A1,n,5,B1,x1,min_eps);A2[0]=2*x1[0],A2[1]=A2[2]=x1[2],A2[3]=2*x1[1];B2[0]=x1[3],B2[1]=x1[4];float *x2=new float[2];//标准化，将一次项消掉，求出center.x和center.y;SVD(A2,2,2,B2,x2,min_eps);rp[0]=x2[0],rp[1]=x2[1];for( i = 0; i < n; i++ ){float px,py;px = arrayx[i]*1.0;py = arrayy[i]*1.0;px -= cx;py -= cy;B3[i] = 1.0;int step=i*3;A3[step] = (px - rp[0]) * (px - rp[0]);A3[step+1] = (py - rp[1]) * (py - rp[1]);A3[step+2] = (px - rp[0]) * (py - rp[1]);}//求出A(x-center.x)^2+B(y-center.y)^2+C(x-center.x)(y-center.y)的最小二乘解SVD(A3,n,3,B3,x1,min_eps);rp[4] = -0.5 * atan2(x1[2], x1[1] - x1[0]);t = sin(-2.0 * rp[4]);if( fabs(t) > fabs(x1[2])*min_eps )t = x1[2]/t;elset = x1[1] - x1[0];rp[2] = fabs(x1[0] + x1[1] - t);if( rp[2] > min_eps )rp[2] = sqrt(2.0 / rp[2]);rp[3] = fabs(x1[0] + x1[1] + t);if( rp[3] > min_eps )rp[3] = sqrt(2.0 / rp[3]);box[0] = (float)rp[0] + cx;box[1]= (float)rp[1] + cy;box[2]= (float)(rp[2]*2);box[3] = (float)(rp[3]*2);if( box[2] > box[3] ){double tmp=box[2];box[2]=box[3];box[3]=tmp;}box[4] = (float)(90 + rp[4]*180/3.1415926);if( box[4] < -180 )box[4] += 360;if( box[4] > 360 )box[4] -= 360;delete []A1;delete []A2;delete []A3;delete []B1;delete []B2;delete []B3;delete []x1;delete []x2;}int LSEllipse::SVD(float *a,int m,int n,float b[],float x[],float esp)
{  float *aa;float *u;float *v;aa=new float[n*m];u=new  float[m*m];v=new  float[n*n];int ka;int  flag;if(m>n){ ka=m+1;}else{ka=n+1;}flag=gmiv(a,m,n,b,x,aa,esp,u,v,ka);delete []aa;delete []u;delete []v;return(flag);
}int LSEllipse::gmiv( float a[],int m,int n,float b[],float x[],float aa[],float eps,float u[],float v[],int ka)
{ int i,j;i=ginv(a,m,n,aa,eps,u,v,ka);if (i<0) return(-1);for (i=0; i<=n-1; i++){ x[i]=0.0;for (j=0; j<=m-1; j++)x[i]=x[i]+aa[i*m+j]*b[j];}return(1);}int LSEllipse::ginv(float a[],int m,int n,float aa[],float eps,float u[],float v[],int ka){ //  int muav(float a[],int m,int n,float u[],float v[],float eps,int ka);int i,j,k,l,t,p,q,f;i=muav(a,m,n,u,v,eps,ka);if (i<0) return(-1);j=n;if (m<n) j=m;j=j-1;k=0;while ((k<=j)&&(a[k*n+k]!=0.0)) k=k+1;k=k-1;for (i=0; i<=n-1; i++)for (j=0; j<=m-1; j++){ t=i*m+j; aa[t]=0.0;for (l=0; l<=k; l++){ f=l*n+i; p=j*m+l; q=l*n+l;aa[t]=aa[t]+v[f]*u[p]/a[q];}}return(1);}int LSEllipse::muav(float a[],int m,int n,float u[],float v[],float eps,int ka){ int i,j,k,l,it,ll,kk,ix,iy,mm,nn,iz,m1,ks;float d,dd,t,sm,sm1,em1,sk,ek,b,c,shh,fg[2],cs[2];float *s,*e,*w;//void ppp();// void sss();void ppp(float a[],float e[],float s[],float v[],int m,int n);void sss(float fg[],float cs[]);s=(float *) malloc(ka*sizeof(float));e=(float *) malloc(ka*sizeof(float));w=(float *) malloc(ka*sizeof(float));it=60; k=n;if (m-1<n) k=m-1;l=m;if (n-2<m) l=n-2;if (l<0) l=0;ll=k;if (l>k) ll=l;if (ll>=1){ for (kk=1; kk<=ll; kk++){ if (kk<=k){ d=0.0;for (i=kk; i<=m; i++){ ix=(i-1)*n+kk-1; d=d+a[ix]*a[ix];}s[kk-1]=(float)sqrt(d);if (s[kk-1]!=0.0){ ix=(kk-1)*n+kk-1;if (a[ix]!=0.0){ s[kk-1]=(float)fabs(s[kk-1]);if (a[ix]<0.0) s[kk-1]=-s[kk-1];}for (i=kk; i<=m; i++){ iy=(i-1)*n+kk-1;a[iy]=a[iy]/s[kk-1];}a[ix]=1.0f+a[ix];}s[kk-1]=-s[kk-1];}if (n>=kk+1){ for (j=kk+1; j<=n; j++){ if ((kk<=k)&&(s[kk-1]!=0.0)){ d=0.0;for (i=kk; i<=m; i++){ ix=(i-1)*n+kk-1;iy=(i-1)*n+j-1;d=d+a[ix]*a[iy];}d=-d/a[(kk-1)*n+kk-1];for (i=kk; i<=m; i++){ ix=(i-1)*n+j-1;iy=(i-1)*n+kk-1;a[ix]=a[ix]+d*a[iy];}}e[j-1]=a[(kk-1)*n+j-1];}}if (kk<=k){ for (i=kk; i<=m; i++){ ix=(i-1)*m+kk-1; iy=(i-1)*n+kk-1;u[ix]=a[iy];}}if (kk<=l){ d=0.0;for (i=kk+1; i<=n; i++)d=d+e[i-1]*e[i-1];e[kk-1]=(float)sqrt(d);if (e[kk-1]!=0.0){ if (e[kk]!=0.0){ e[kk-1]=(float)fabs(e[kk-1]);if (e[kk]<0.0) e[kk-1]=-e[kk-1];}for (i=kk+1; i<=n; i++)e[i-1]=e[i-1]/e[kk-1];e[kk]=1.0f+e[kk];}e[kk-1]=-e[kk-1];if ((kk+1<=m)&&(e[kk-1]!=0.0)){ for (i=kk+1; i<=m; i++) w[i-1]=0.0;for (j=kk+1; j<=n; j++)for (i=kk+1; i<=m; i++)w[i-1]=w[i-1]+e[j-1]*a[(i-1)*n+j-1];for (j=kk+1; j<=n; j++)for (i=kk+1; i<=m; i++){ ix=(i-1)*n+j-1;a[ix]=a[ix]-w[i-1]*e[j-1]/e[kk];}}for (i=kk+1; i<=n; i++)v[(i-1)*n+kk-1]=e[i-1];}}}mm=n;if (m+1<n) mm=m+1;if (k<n) s[k]=a[k*n+k];if (m<mm) s[mm-1]=0.0;if (l+1<mm) e[l]=a[l*n+mm-1];e[mm-1]=0.0;nn=m;if (m>n) nn=n;if (nn>=k+1){ for (j=k+1; j<=nn; j++){ for (i=1; i<=m; i++)u[(i-1)*m+j-1]=0.0;u[(j-1)*m+j-1]=1.0;}}if (k>=1){ for (ll=1; ll<=k; ll++){ kk=k-ll+1; iz=(kk-1)*m+kk-1;if (s[kk-1]!=0.0){ if (nn>=kk+1)for (j=kk+1; j<=nn; j++){ d=0.0;for (i=kk; i<=m; i++){ ix=(i-1)*m+kk-1;iy=(i-1)*m+j-1;d=d+u[ix]*u[iy]/u[iz];}d=-d;for (i=kk; i<=m; i++){ ix=(i-1)*m+j-1;iy=(i-1)*m+kk-1;u[ix]=u[ix]+d*u[iy];}}for (i=kk; i<=m; i++){ ix=(i-1)*m+kk-1; u[ix]=-u[ix];}u[iz]=1.0f+u[iz];if (kk-1>=1)for (i=1; i<=kk-1; i++)u[(i-1)*m+kk-1]=0.0;}else{ for (i=1; i<=m; i++)u[(i-1)*m+kk-1]=0.0;u[(kk-1)*m+kk-1]=1.0;}}}for (ll=1; ll<=n; ll++){ kk=n-ll+1; iz=kk*n+kk-1;if ((kk<=l)&&(e[kk-1]!=0.0)){ for (j=kk+1; j<=n; j++){ d=0.0;for (i=kk+1; i<=n; i++){ ix=(i-1)*n+kk-1; iy=(i-1)*n+j-1;d=d+v[ix]*v[iy]/v[iz];}d=-d;for (i=kk+1; i<=n; i++){ ix=(i-1)*n+j-1; iy=(i-1)*n+kk-1;v[ix]=v[ix]+d*v[iy];}}}for (i=1; i<=n; i++)v[(i-1)*n+kk-1]=0.0;v[iz-n]=1.0;}for (i=1; i<=m; i++)for (j=1; j<=n; j++)a[(i-1)*n+j-1]=0.0;m1=mm; it=60;while (1==1){ if (mm==0){ ppp(a,e,s,v,m,n);free(s); free(e); free(w); return(1);}if (it==0){ ppp(a,e,s,v,m,n);free(s); free(e); free(w); return(-1);}kk=mm-1;while ((kk!=0)&&(fabs(e[kk-1])!=0.0)){ d=(float)(fabs(s[kk-1])+fabs(s[kk]));dd=(float)fabs(e[kk-1]);if (dd>eps*d) kk=kk-1;else e[kk-1]=0.0;}if (kk==mm-1){ kk=kk+1;if (s[kk-1]<0.0){ s[kk-1]=-s[kk-1];for (i=1; i<=n; i++){ ix=(i-1)*n+kk-1; v[ix]=-v[ix];}}while ((kk!=m1)&&(s[kk-1]<s[kk])){ d=s[kk-1]; s[kk-1]=s[kk]; s[kk]=d;if (kk<n)for (i=1; i<=n; i++){ ix=(i-1)*n+kk-1; iy=(i-1)*n+kk;d=v[ix]; v[ix]=v[iy]; v[iy]=d;}if (kk<m)for (i=1; i<=m; i++){ ix=(i-1)*m+kk-1; iy=(i-1)*m+kk;d=u[ix]; u[ix]=u[iy]; u[iy]=d;}kk=kk+1;}it=60;mm=mm-1;}else{ ks=mm;while ((ks>kk)&&(fabs(s[ks-1])!=0.0)){ d=0.0;if (ks!=mm) d=d+(float)fabs(e[ks-1]);if (ks!=kk+1) d=d+(float)fabs(e[ks-2]);dd=(float)fabs(s[ks-1]);if (dd>eps*d) ks=ks-1;else s[ks-1]=0.0;}if (ks==kk){ kk=kk+1;d=(float)fabs(s[mm-1]);t=(float)fabs(s[mm-2]);if (t>d) d=t;t=(float)fabs(e[mm-2]);if (t>d) d=t;t=(float)fabs(s[kk-1]);if (t>d) d=t;t=(float)fabs(e[kk-1]);if (t>d) d=t;sm=s[mm-1]/d; sm1=s[mm-2]/d;em1=e[mm-2]/d;sk=s[kk-1]/d; ek=e[kk-1]/d;b=((sm1+sm)*(sm1-sm)+em1*em1)/2.0f;c=sm*em1; c=c*c; shh=0.0;if ((b!=0.0)||(c!=0.0)){ shh=(float)sqrt(b*b+c);if (b<0.0) shh=-shh;shh=c/(b+shh);}fg[0]=(sk+sm)*(sk-sm)-shh;fg[1]=sk*ek;for (i=kk; i<=mm-1; i++){ sss(fg,cs);if (i!=kk) e[i-2]=fg[0];fg[0]=cs[0]*s[i-1]+cs[1]*e[i-1];e[i-1]=cs[0]*e[i-1]-cs[1]*s[i-1];fg[1]=cs[1]*s[i];s[i]=cs[0]*s[i];if ((cs[0]!=1.0)||(cs[1]!=0.0))for (j=1; j<=n; j++){ ix=(j-1)*n+i-1;iy=(j-1)*n+i;d=cs[0]*v[ix]+cs[1]*v[iy];v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];v[ix]=d;}sss(fg,cs);s[i-1]=fg[0];fg[0]=cs[0]*e[i-1]+cs[1]*s[i];s[i]=-cs[1]*e[i-1]+cs[0]*s[i];fg[1]=cs[1]*e[i];e[i]=cs[0]*e[i];if (i<m)if ((cs[0]!=1.0)||(cs[1]!=0.0))for (j=1; j<=m; j++){ ix=(j-1)*m+i-1;iy=(j-1)*m+i;d=cs[0]*u[ix]+cs[1]*u[iy];u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];u[ix]=d;}}e[mm-2]=fg[0];it=it-1;}else{ if (ks==mm){ kk=kk+1;fg[1]=e[mm-2]; e[mm-2]=0.0;for (ll=kk; ll<=mm-1; ll++){ i=mm+kk-ll-1;fg[0]=s[i-1];sss(fg,cs);s[i-1]=fg[0];if (i!=kk){ fg[1]=-cs[1]*e[i-2];e[i-2]=cs[0]*e[i-2];}if ((cs[0]!=1.0)||(cs[1]!=0.0))for (j=1; j<=n; j++){ ix=(j-1)*n+i-1;iy=(j-1)*n+mm-1;d=cs[0]*v[ix]+cs[1]*v[iy];v[iy]=-cs[1]*v[ix]+cs[0]*v[iy];v[ix]=d;}}}else{ kk=ks+1;fg[1]=e[kk-2];e[kk-2]=0.0;for (i=kk; i<=mm; i++){ fg[0]=s[i-1];sss(fg,cs);s[i-1]=fg[0];fg[1]=-cs[1]*e[i-1];e[i-1]=cs[0]*e[i-1];if ((cs[0]!=1.0)||(cs[1]!=0.0))for (j=1; j<=m; j++){ ix=(j-1)*m+i-1;iy=(j-1)*m+kk-2;d=cs[0]*u[ix]+cs[1]*u[iy];u[iy]=-cs[1]*u[ix]+cs[0]*u[iy];u[ix]=d;}}}}}}free(s);free(e);free(w); return(1);}void ppp(float a[],float e[],float s[],float v[],int m,int n)
{ int i,j,p,q;float d;if (m>=n) i=n;else i=m;for (j=1; j<=i-1; j++){ a[(j-1)*n+j-1]=s[j-1];a[(j-1)*n+j]=e[j-1];}a[(i-1)*n+i-1]=s[i-1];if (m<n) a[(i-1)*n+i]=e[i-1];for (i=1; i<=n-1; i++)for (j=i+1; j<=n; j++){ p=(i-1)*n+j-1; q=(j-1)*n+i-1;d=v[p]; v[p]=v[q]; v[q]=d;}return;}void sss(float fg[],float cs[]){ float r,d;if ((fabs(fg[0])+fabs(fg[1]))==0.0){ cs[0]=1.0; cs[1]=0.0; d=0.0;}else { d=(float)sqrt(fg[0]*fg[0]+fg[1]*fg[1]);if (fabs(fg[0])>fabs(fg[1])){ d=(float)fabs(d);if (fg[0]<0.0) d=-d;}if (fabs(fg[1])>=fabs(fg[0])){ d=(float)fabs(d);if (fg[1]<0.0) d=-d;}cs[0]=fg[0]/d; cs[1]=fg[1]/d;}r=1.0;if (fabs(fg[0])>fabs(fg[1])) r=cs[1];elseif (cs[0]!=0.0) r=1.0f/cs[0];fg[0]=d; fg[1]=r;return;}

参考:

线性，非线性多项式：

http://blog.csdn.net/qll125596718/article/details/8248249

http://blog.csdn.net/poxiaozhuimeng/article/details/41117947

http://blog.csdn.net/ouyangying123/article/details/53996403

http://blog.csdn.net/jairuschan/article/details/7517773/

http://blog.csdn.net/zang141588761/article/details/50523036

http://www.cnblogs.com/gnuhpc/archive/2012/12/09/2809997.html

http://download.csdn.net/download/biaobiao11/9755119

圆拟合：

http://blog.sina.com.cn/s/blog_b27f71160101gxun.html

http://www.cnblogs.com/dotLive/archive/2007/04/06/524633.html

http://blog.csdn.net/andylao62/article/details/24522365

http://blog.csdn.net/liyuanbhu/article/details/50889951

http://blog.csdn.net/liyuanbhu/article/details/50890587

椭圆拟合：

http://doc.okbase.net/u013708970/archive/121532.html