matlib.cpp

基于线性规划的回归支持向量机源程序

                        if ((cs[0]!=1.0)||(cs[1]!=0.0))
							for (j=1; j<=m; j++){
 								d=cs[0]*u[j-1][i-1]+cs[1]*u[j-1][kk-2];
                                u[j-1][kk-2]=-cs[1]*u[j-1][i-1]+cs[0]*u[j-1][kk-2];
								u[j-1][i-1]=d;
                            }
					}
				}
			}
          }
      }
}
matrix pinv(matrix& a,double eps)
{
    int m=a.mrow(),n=a.mcol();
    matrix aa(n,m);
    double data=0;
    int k;
    if(m==1)
    {
        data=dot(a[0],a[0]);
        for(k=0;k<n;k++)
        aa[k][0]=a[0][k];
        if(data>0.000000001)
        aa/=data;
        return aa;
    }
    else
    if(n==1)
    {
        for(k=0;k<m;k++)
        aa[0][k]=a[k][0];
        data=dot(aa[0],aa[0]);
        if(data>0.000000001)
        aa/=data;
        return aa;
    }
    matrix u(m,m),v(n,n),x;
	x=a.copy();
    int i,j=MIN(m,n)-1,l;
    svd(x,u,v,eps);
	k=0;
    while ((k<=j)&&(fabs(x[k][k])>0.000000001)) k++;
	k--;
    for (i=0; i<n; i++)
		for (j=0; j<m; j++)
			for (l=0; l<=k; l++)
				aa[i][j]+=v[l][i]*u[j][l]/x[l][l];
    return aa;
}
int dcinv(matrix& a)
{
	int n=mrow(a);
	IntVector is(n);
	IntVector js(n);
    double d,p;
	int i,j,k;
    for (k=0; k<=n-1; k++)
    {
		d=0.0;
        for (i=k; i<=n-1; i++)
        for (j=k; j<=n-1; j++)
        {
			p=fabs(a(i,j));
            if (p>d) { d=p; is[k]=i; js[k]=j;}
        }
        if (d+1.0==1.0)
        {
			return 0;
        }
        if (is[k]!=k)
			for (j=0; j<=n-1; j++)
            {
				p=a(k,j);
				a(k,j)=a(is[k],j);
				a(is[k],j)=p;
            }
        if (js[k]!=k)
			for (i=0; i<=n-1; i++)
            {
				p=a(i,k); a(i,k)=a(i,js[k]); a(i,js[k])=p;
            }
        p=1.0/a(k,k);
		a(k,k)=p;
        for (j=0; j<=n-1; j++)
			if (j!=k) a(k,j)*=p;
        for (i=0; i<=n-1; i++)
          if (i!=k)
            for (j=0; j<=n-1; j++)
	            if (j!=k) a(i,j)-=(a(i,k)*a(k,j));
        for (i=0; i<=n-1; i++)
	        if (i!=k) a(i,k)*=(-p);
      }
    for (k=n-1; k>=0; k--)
    {
		if (js[k]!=k)
          for (j=0; j<=n-1; j++)
            {
				p=a(k,j); a(k,j)=a(js[k],j); a(js[k],j)=p;
            }
        if (is[k]!=k)
          for (i=0; i<=n-1; i++)
            {
              p=a(i,k); a(i,k)=a(i,is[k]); a(i,is[k])=p;
            }
      }
    return 1;
}
// max f=c'*x subject to: A*x<=b,b>=0,x>=0;
int lp_simplex(matrix& a, vector& b,vector& c, vector& x)
{
	int i,k,j,l;
	int m=mrow(a);
	int n=mcol(a);
	x.initial(n);
    double s,z,dd,y;
    IntVector js;
	js.initial(m);
    matrix p(m,m);
    matrix d(m,n);
	x.initial(n);
    for (i=0; i<m; i++) js[i]=n-m+i;
    s=0.0;
    while (1==1)
    { 
		for (i=0; i<m; i++)
			for (j=0; j<m; j++)
				p(i,j)=a(i,js[j]);
        l=dcinv(p);
        if (l==0)
        { 
			x[n-m]=s;
            return l;
        }
        d=p*a;
        for (i=0; i<=n-1; i++) x[i]=0.0;
        for (i=0; i<=m-1; i++)
        {
			s=0.0;
            for (j=0; j<=m-1; j++)
			    s=s+p(i,j)*b[j];
            x[js[i]]=s;
        }
        k=-1; dd=1.0e-35;
        for (j=0; j<=n-1; j++)
        {
			z=0.0;
            for (i=0; i<=m-1; i++)
            z=z+c[js[i]]*d(i,j);
            z=z-c[j];
            if (z>dd) { dd=z; k=j;}
          }
        if (k==-1)
        {
			s=0.0;
            for (j=0; j<=n-m-1; j++)
              s=s+c[j]*x[j];
			x[n-m]=s;
            return 1;
        }
        j=-1;
        dd=DOUBLE_MAX;//1.0e+20;
        for (i=0; i<=m-1; i++)
        if (d(i,k)>=1.0e-20)
            { y=x[js[i]]/d(i,k);
              if (y<dd) { dd=y; j=i;}
            }
        if (j==-1) 
		{
			x[n-m]=s;
			return 0;
        }
        js[j]=k;
	}
	return 1;
}
/*标准单纯形法，直接求逆矩阵
int lp_simplex(matrix& a, vector& b,vector& c, IntVector& js, vector& x)
{
	int i,k,j,l;
	int m=mrow(a);
	int n=mcol(a);
	x.initial(n);
    double s,z,dd,y;
    matrix p(m,m);
    matrix d(m,n);
	x.initial(n+1);
	IntVector is(n-m);
	k=0;
	for(i=0;i<n;i++)
	{
		int f=0;
		for(j=0;j<m;j++)
		{
			if(i==js[j])
			{
				f=1;
				break;
			}
		}
		if(f)
			continue;
		is[k]=i;
		k++;
	}
    s=0.0;
    while (1==1)
    { 
		for (i=0; i<m; i++)
			for (j=0; j<m; j++)
				p(i,j)=a(i,js[j]);
        l=dcinv(p);
        if (l==0)
        { 
			x[n]=s;
            return l;
        }
        d=p*a;
        x=0.0;
        for (i=0; i<=m-1; i++)
        {
			s=0.0;
            for (j=0; j<=m-1; j++)
			    s=s+p(i,j)*b[j];
            x[js[i]]=s;
        }
        k=-1; dd=1.0e-13;
        for (j=0; j<=n-m-1; j++)
        {
			z=0.0;
            for (i=0; i<=m-1; i++)
            z=z+c[js[i]]*d(i,is[j]);
            z=z-c[is[j]];
            if (z>dd) 
			{ dd=z; k=j;}
        }
        if (k==-1)
        {
			s=0.0;
            for (j=0; j<=n-1; j++)
              s=s+c[j]*x[j];
			x[n]=s;
            return 1;
        }
        j=-1;
        dd=DOUBLE_MAX;
        for (i=0; i<=m-1; i++)
        if (d(i,is[k])>=1.0e-12)
            { y=x[js[i]]/d(i,is[k]);
              if (y<dd) { dd=y; j=i;}
            }
        if (j==-1) 
		{
			x[n]=s;
			return 0;
        }
        i=js[j];
		js[j]=is[k];
		is[k]=i;
//cout<<"d("<<j<<","<<k<<")="<<d(j,js[j])<<endl;
//getch();
	}
	return 1;
}
*/
//标准单纯形法，采用一列置换求逆算法计算逆矩阵
int lp_simplex(matrix& a, vector& b,vector& c, IntVector& js, vector& x)
{
	int i,k,j,l;
	int m=mrow(a);
	int n=mcol(a);
	x.initial(n);
    double s,z,dd,y;
    matrix p(m,m);
    matrix d(m,n);
	x.initial(n+1);
	IntVector is(n-m);
	k=0;
	for(i=0;i<n;i++)
	{
		int f=0;
		for(j=0;j<m;j++)
		{
			if(i==js[j])
			{
				f=1;
				break;
			}
		}
		if(f)
			continue;
		is[k]=i;
		k++;
	}
    s=0.0;
	for (i=0; i<m; i++)
		for (j=0; j<m; j++)
			p(i,j)=a(i,js[j]);
    l=dcinv(p);
    while (1)
    { 
        if (l==0)
        { 
			x[n]=s;
            return l;
        }
        d=p*a;
        d=p*a;
        x=0.0;
        for (i=0; i<=m-1; i++)
        {
			s=0.0;
            for (j=0; j<=m-1; j++)
			    s=s+p(i,j)*b[j];
            x[js[i]]=s;
        }
        k=-1; dd=1.0e-13;
        for (j=0; j<=n-m-1; j++)
        {
			z=0.0;
            for (i=0; i<=m-1; i++)
            z=z+c[js[i]]*d(i,is[j]);
            z=z-c[is[j]];
            if (z>dd) 
			{ dd=z; k=j;}
        }
        if (k==-1)
        {
			s=0.0;
            for (j=0; j<=n-1; j++)
              s=s+c[j]*x[j];
			x[n]=s;
            return 1;
        }
        j=-1;
        dd=DOUBLE_MAX;
        for (i=0; i<=m-1; i++)
        if (d(i,is[k])>=1.0e-12)
            { y=x[js[i]]/d(i,is[k]);
              if (y<dd) { dd=y; j=i;}
            }
        if (j==-1) 
		{
			x[n]=s;
			return 0;
        }
		double data=d(j,is[k]);
		if(fabs(data)<=0.00000000001)
		{
			x[n]=s;
			return -1;
		}
		d(j,is[k])-=1;
		for(i=0;i<m;i++)
		d(i,js[j])=p(j,i);
		for(i=0;i<m;i++)
		for(int t=0;t<m;t++)
		p(i,t)-=d(i,is[k])*d(t,js[j])/data;

        i=js[j];
		js[j]=is[k];
		is[k]=i;
//cout<<"d("<<j<<","<<k<<")="<<d(j,js[j])<<endl;
//getch();
	}
	return 1;
}

int lp_affineinner(matrix& a, vector& b,vector& c, vector& x)
{
	int i,k,j,l;
	int fail=0;
	int m=mrow(a);
	int n=mcol(a);
	double e1=0.1;
	double e2=1.0e-10;
	double e3=1.0e-10;
	int f1,f2,f3;
	x.initial(n);
	vector w(m),s(n);
	w=0;x=1;s=1;
	double bn1=0,cn1=0;
	for(i=0;i<b();i++)
	{
		double data=fabs(b[i]);
		if(bn1<data)
			bn1=data;
	}
	bn1+=1;
	for(i=0;i<c();i++)
	{
		double data=fabs(c[i]);
		if(cn1<data)
			cn1=data;
	}
	cn1+=1;
	vector u,t,v,p,d2;
	vector dw,ds,dx;
	matrix D(m,m);
	k=1;
	while(1)
	{
		double mu;
		mu=dot(x,s)/n;
		t=b-a*x;
		u=c-w*a-s;
		v=mu-x*s;
		p=v/x;
		d2=x/s;
//cout<<"t="<<t<<endl;
//cout<<"u="<<u<<endl;
//cout<<"v="<<v<<endl;
//cout<<"p="<<p<<endl;
		if(mu<e1)
			f1=1;
		else
			f1=0;
		double data1=0,data2=0;
		for(i=0;i<t();i++)
		{
			double data=fabs(t[i]);
			if(data1<data)
				data1=data;
		}
		for(i=0;i<u();i++)
		{
			double data=fabs(u[i]);
			if(data2<data)
				data2=data;
		}
//cout<<"x="<<x<<endl;
//cout<<"w="<<w<<endl;
//cout<<"s="<<s<<endl;
cout<<"mu="<<mu<<", tmax="<<data1<<", umax="<<data2<<endl;
		data1/=bn1;
		data2/=cn1;
		if(data1<=e2)
			f2=1;
		else
			f2=0;
		if(data2<=e3)
			f3=1;
		else
		f3=0;
		if(f1&&f2&&f3)
			break;
		vector up;
		up=u-p;
		up*=d2;
		v=a*up+t;
		for(i=0;i<m;i++)
			for(j=i;j<m;j++)
			{
				double dd=0;
				for(l=0;l<n;l++)
				{
					double data=a[i][l]*d2[l]*a[j][l];
					dd+=data;
				}
				D(i,j)=dd;
				if(i!=j)
					D(j,i)=dd;
			}
//cout<<"D=\n"<<D;
//cout<<"v="<<v<<endl;
		SolveSysEqa(D,v,dw);
		ds=u-dw*a;
		dx=p-ds;
		dx*=d2;
//cout<<"dw=inv(D)*v="<<dw<<endl;
//cout<<"ds="<<ds<<endl;
//cout<<"dx="<<dx<<endl;
		if(f2)
		{
			data1=dx.vmin();
			data2=dx.vmax();
			if((data1>=-1.0e-20)&&(data2>=e2))
				f1=1;
			else
				f1=0;
			data1=dot(c,dx);
			if( (data1<1.0e-20) && f1)
				fail=1;
		}
		if(f3)
		{
			data1=ds.vmin();
			data2=ds.vmax();
			if((data1>=-1.0e-20)&&(data2>=e3))
				f2=1;
			else
				f2=0;
			data1=dot(b,dw);
			if( (data1>e3) && f2)
				fail=-1;
		}
		if(fail)
			break;
		data1=1;
		data2=1;
		double alpha=0.99;
		for(i=0;i<n;i++)
		{
			double data=-dx[i]/(alpha*x[i]);
			if(data1<data)
				data1=data;
			data=-ds[i]/(alpha*s[i]);
			if(data2<data)
				data2=data;
		}
cout<<"px="<<1/data1<<",  pw="<<1/data2<<endl;
if(data1>data2)
data1=data2;
		dx/=data1;
		dw/=data2;
		ds/=data2;
		x+=dx;
		w+=dw;
		s+=ds;
k++;
cout<<"\n k="<<k<<endl;
getch();
	}
	return fail;
}
int dicll(matrix& a)
{
	int n=a.mrow();
    double d;
    if ((a(0,0)+1.0==1.0)||(a(0,0)<0.0))
      { printf("fail\n"); return(-2);}
    a(0,0)=sqrt(a(0,0));
    d=a(0,0);
    for(int i=1; i<=n-1; i++)
    {
		a(i,0)=a(i,0)/d;
	}
    for (int j=1; j<=n-1; j++)
    {
		double data=a(j,j);
        for(int k=0; k<=j-1; k++)
        {
			data=data-a(j,k)*a(j,k);
		}
        if((data+1.0==1.0)||(data<0.0))
          { printf("fail\n"); return(-2);}
        a(j,j)=sqrt(data);
        d*=a(j,j);
        for(i=j+1; i<=n-1; i++)
        { 
            for (k=0; k<=j-1; k++)
            a(i,j)-=a(i,k)*a(j,k);
            a(i,j)/=a(j,j);
        }
      }
    for (i=0; i<=n-2; i++)
    for (j=i+1; j<=n-1; j++)
      a(i,j)=0.0;
    return 1;
}

//mrow(a)=vlen(v),a*x=v,vlen(x)=mcol(a)

int SolveSysEqa(matrix& a,vector& v,vector& x)
{
    int m=mrow(a);
    int n=mcol(a);
    if((m!=n)||(n!=v()))
    return 0;
	int f;
	f=dicll(a);
	if(f<=0)
		return 0;
	x.initial(m);
	vector u(m);
	u[0]=v[0]/a(0,0);
	for(int i=1;i<m;i++)
	{
		double data=v[i];
		for(int j=0;j<i;j++)
			data-=a(i,j)*u[j];
		u[i]=data/a(i,i);
	}
  	x[m-1]=u[m-1]/a(m-1,m-1);
	for(i=m-2;i>=0;i--)
	{
		double data=u[i];
		for(int j=i+1;j<m;j++)
			data-=a(j,i)*x[j];
		x[i]=data/a(i,i);
	}
  return 1;
}
int lu(matrix& a,IntVector& js)
{
/*
	int mr=mrow(a);
	int mc=mcol(a);
	js.initial(mr);
	IntVector f(mr);
	for(int i=0;i<mr;i++)
		js[i]=i;
	f=1;
	matrix p(mr,mr);
	for(i=0;i<mr;i++)
	{
		double d=fabs(a(i,i));
		int s=i;
		for(int j=i+1;j<mr;j++)
		{
			double data=fabs(a(i,j));
			if(d<data)
			{
				d=data;
				s=j;
			}
		}
		if(s>i)
		{
			js[i]=s;
			js[s]=i;
		}
		for(int j=0;j<i;j++)
		{
			double data=a(js[i],j);
			for(int k=0;k<j;k++)
				data-=p(i,k)*p(k,j);
			data/=p(j,j);
		}

	}
*/
return 0;
}