最小二乘法记录

最小二乘法记录

模板：1668. Death Star 2 @ Timus Online Judge

最小二乘法是计算一个一个列向量\(x\)满足以下式子最小：

\(|A\times x-b|\)，其中\(A\)是\(n\)行\(m\)列矩阵，\(x\)是\(m\)行向量，\(b\)是\(n\)行向量。

也就是上面的模板题。

写成清楚的形式：

\[\min{\sqrt{\sum_{i=1}^n[(\sum_{j=1}^m a_{i,j}x_j)-b_i]^2}} \]

求\(x\)。

方法：

首先把\(\sqrt {...}\)去掉，然后需要求\(\sum_{i=1}^n[(\sum_{j=1}^m a_{i,j}x_j)-b_i]^2\)取到极值的\(x\)（通常来讲，这个式子有下届无上界）。而取到极值是，导数为0(可以理解成一个切面，如果是顶峰就是平面)，所以我们对于每一个\(x_k\)求出偏导，记作\(g_k(x)\)。

\[g_k(x)=\sum_{i=1}^n2a_{i,k}\sum a_{i,j}x_j-2b_ia_{i,k}=\sum _{i=1}^n2a_{i,k}(\sum_{j=1}^m a_{i,j}x_j-b_i) \]

可以写成很简单的矩阵形式：

\[g_k(x)=2[A^T(Ax-b)]_k \]

也就是说，需要找到一个\(x\),使得\(\forall_{1\leq k\leq m}g_k(x)=0\)。

直接高斯消元即可。通常来讲，任意一组解就是最优解，但是我也不会证明。

例题代码 ural 1668:

/*
{
######################
#       Author       #
#        Gary        #
#        2021        #
######################
*/
#include<bits/stdc++.h>
#define rb(a,b,c) for(int a=b;a<=c;++a)
#define rl(a,b,c) for(int a=b;a>=c;--a)
#define LL long long
#define IT iterator
#define PB push_back
#define II(a,b) make_pair(a,b)
#define FIR first
#define SEC second
#define FREO freopen("check.out","w",stdout)
#define rep(a,b) for(int a=0;a<b;++a)
#define SRAND mt19937 rng(chrono::steady_clock::now().time_since_epoch().count())
#define random(a) rng()%a
#define ALL(a) a.begin(),a.end()
#define POB pop_back
#define ff fflush(stdout)
#define fastio ios::sync_with_stdio(false)
#define check_min(a,b) a=min(a,b)
#define check_max(a,b) a=max(a,b)
using namespace std;
//inline int read(){
//    int x=0;
//    char ch=getchar();
//    while(ch<'0'||ch>'9'){
//        ch=getchar();
//    }
//    while(ch>='0'&&ch<='9'){
//        x=(x<<1)+(x<<3)+(ch^48);
//        ch=getchar();
//    }
//    return x;
//}
const int INF=0x3f3f3f3f;
typedef pair<int,int> mp;
/*}
*/
const double eps=1e-8;
namespace GAUSS{
	void gauss(vector<vector<double> >& mat,int n,int m){
		rb(i,1,m){
			int x=-1;
			rb(j,i,n){
				if(abs(mat[j][i])>eps){
					x=j;
					break;
				}
			}
			if(x==-1) continue;
			swap(mat[x],mat[i]);
			double tmp=1/mat[i][i];
			rb(j,1,m+1) mat[i][j]*=tmp;
			rb(j,1,n){
				if(j!=i&&abs(mat[j][i])>eps){
					double val=mat[j][i];			
					rb(k,1,m+1){
						mat[j][k]-=mat[i][k]*val;
					}
				}
			}
		}
	}
}
const int MAXN=101;
int n,m,a[MAXN][MAXN],b[MAXN];
bool fre[MAXN];
vector<vector<double> > mat;
vector<vector<double> > mat2;
double rest[MAXN];
int main(){
	scanf("%d%d",&n,&m);
	rb(i,1,m) rb(j,1,n) scanf("%d",&a[j][i]);
	rb(i,1,n) scanf("%d",&b[i]);
	mat.resize(m+1);
	rb(i,1,m) mat[i].resize(m+2);
	mat2=mat;
	rb(k,1,m) rb(i,1,n){
		rb(j,1,m) mat[k][j]+=2.0*a[i][k]*a[i][j];
		mat[k][m+1]+=2.0*b[i]*a[i][k];
	}
	GAUSS:: gauss(mat,m,m);
	rb(i,1,m){
		if(abs(mat[i][i])<eps){
			mat[i][i]=1.0;
		}
		else{
			mat[i][i]=0.0;	
			rb(j,1,m+1) mat[i][j]*=-1;
		}
	}
	rb(k,1,m) rb(i,1,m){
		rb(j,1,m) mat2[k][j]+=2.0*mat[i][k]*mat[i][j];
		mat2[k][m+1]+=2.0*mat[i][m+1]*mat[i][k]; 
	}
	GAUSS:: gauss(mat2,m,m);
	rb(i,1,m) if(abs(mat2[i][i])<eps) fre[i]=true;
	rb(i,1,m) if(fre[i]) rest[i]=114.514;
	rb(i,1,m) if(!fre[i]){
		rb(j,1,m) if(fre[j]) mat2[i][m+1]-=mat2[i][j]*rest[j];
		rest[i]=mat2[i][m+1];
	}
	rb(i,1,m){
		double ans=0;
		rb(j,1,m){
			ans+=mat[i][j]*rest[j];
		}
		ans-=mat[i][m+1];
		printf("%.5f ",ans);
	}
	puts("");
	return 0;
}