matlib.cpp

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

    return y;
    y.initial(x());
    for(int i=0;i<x();i++)
    {
        double d=p[0];
        double data=x[i];
        double d1=data;
        for(int j=1;j<m;j++)
        {
            d+=p[j]*data;
            data*=d1;
        }
        y[i]=d;
    }
    return y;
}

/*
double GetLypExp(matrix& a)
{
    int mr=mrow(a);
    int mc=mcol(a);
    double dmax=max(maxofrow(a));
    double dmin=min(minofrow(a));
    double radio=(dmax-dmin)*0.05;
    radio=radio*radio;
    vector v,u;
    int x1=0,y1=1;
    double dist1=0;
    do
    {
        y1=x1+1;
        double d=-1.0;
        for(int i=x1+1;i<mr;i++)
        {
            double data;
            v=a[x1]-a[i];
            data=dot(v,v);
            if(i==(x1+1))
            d=data;
            else
            if(d>data)
            {
                d=data;
                y1=i;
            }
        }
        if(d<radio)
        break;
        if(x1<(mr-2))
        x1++;
        else
        return 1000000.0;
    }
    while(1);
    dist1=d;
    double dist2;
    int x2,y2;
    x2=x1;
    y2=y1;
    dist2=dist1;
    int N=0;
    double lypexp=0;
    while(x2<(mr-1))
    {
        x2++;
        y2++;
        if(y2>=(mr-1))
        {
            v=a[x2]-a[y2];
            dist2=dot(v,v);
            dmin=(log(dist2-dist1))/(x2-x1);
            lypexp=lypexp*N/(N+1)+dmin/(N+1);
            return lypexp;
        }
        v=a[x2]-a[y2];
        dist2=dot(v,v);
        if(dist2>=radio)
        {
            dmin=(log(dist2-dist1))/(x2-x1);
            lypexp=lypexp*N/(N+1)+dmin/(N+1);
            N++;
            x1=x2;
            y1=-1;
            for(i=0;i<(mr-1);i++)
            {
                if(i==x1)
                continue;
                v=a[x1]-a[i];
                dist1=dot(v,v);
                u=a[x2]-a[y1];
                dmin=dot(u,v);
                if(y1<0)
                {
                    y1=i;
                    dist2=dist1;
                    dmax=dmin;
                    continue;
                }
                if( (dist2<radio)&&(dist1<radio) )
                {
                    if(dmin<dmax)
                    {
                        dmax=dmin;
                        dist2=dist1;
                        y1=i;
                    }
                }
                else
                {
                    if(dist1<radio)
                    {
                        dmax=dmin;
                        dist2=dist1;
                        y1=i;
                    }
                }
            }
            if(dist2<radio)
            {
                dist1=dist2;
                y2=y1;
                x2=x1;
            }
            else
            return lypexp;
        }
    }
    return lypexp;
}
*/
int SolveNewtonDir(matrix& a,vector& b)
{
    vector v;
    matrix L;
    v=copy(b);
    int mr,mc;
    mr=a.mrow();
    mc=a.mcol();
    if(mr!=mc)
    return 0;
    if(mr!=b())
    return 0;
/*
    d1=0;d2=0;
    for(int i=0;i<mr;i++)
    if(d1<fabs(a(i,i)))
    d1=fabs(a(i,i));
    for(i=0;i<mr-1;i++)
    for(int j=i+1;j<mc;j++)
    if(d2<fabs(a(i,j)))
    d2=fabs(a(i,j));
    beta=MAX(d1,d2/mr);
    beta=sqrt(beta);
    delta=0.000001;
    L.initial(mr,mr);
    L=0.0;
    for(i=0;i<mr;i++)
    {
        d1=a(i,i);
        for(int j=0;j<i;j++)
        {
            d2=L(i,j);
            d1-=d2*d2;
        }
        d2=sqrt(fabs(d1));
        d1=MAX(delta,d2);
        L(i,i)=d1;
        for(j=i+1;j<mr;j++)
        {
            d2=a(j,i);
            for(int k=0;k<i;k++)
            d2-=L(j,k)*L(i,k);
            L(j,i)=d2/d1;
        }

        d1=0;
        for(j=i+1;j<mr;j++)
        {
            if(d1<fabs(L(j,i)) )
            d1=fabs(L(j,i));
        }
        if(d1>=beta)
        {
            L(i,i)=L(i,i)*d1/beta;
            for(int j=i+1;j<mr;j++)
            L(j,i)=L(j,i)*beta/d1;
        }
    }
    v.initial(mr);
//    cout<<L<<endl;
    for(i=0;i<mr;i++)
    {
        d1=b[i]/L(i,i);
        for(int j=0;j<i;j++)
        d1-=L(i,j)*v[j]/L(i,i);
        v[i]=d1;
    }
    for(i=mr-1;i>=0;i--)
    {
        d1=v[i]/L(i,i);
        for(int j=i+1;j<mr;j++)
        d1-=L(j,i)*v[j]/L(i,i);
        b[i]=d1;
    }


*/
    matrix p;
    SysEig(a,p,v);
    vector u,d;
    u=b*p;
    d.initial(mr);
    for(int i=0;i<mr;i++)
    if(v[i]<0.0)
    {
        if(v[i]<-0.00001)
        d[i]=u[i]/v[i];
        else
        d[i]=-u[i]/0.00001;
    }
    else
    {
        if(v[i]>0.00001)
        d[i]=-u[i]/v[i];
        else
        d[i]=-u[i]/0.00001;
    }
    b=p*d;
    b*=-1.0;

    return 1;
}

int SolveNewtonDir2(matrix& a,vector& b)
{
    vector v;
//    matrix L;
    v=copy(b);
    int mr,mc;
    mr=a.mrow();
    mc=a.mcol();
    if(mr!=mc)
    return 0;
    if(mr!=b())
    return 0;
    matrix p;
    SysMatrixCanonical(a,p,v);
    vector u,d;
    u=b*p;
    d.initial(mr);
    for(int i=0;i<mr;i++)
    if(v[i]<0.0)
    {
        if(v[i]<-0.00001)
        d[i]=u[i]/v[i];
        else
//        d[i]=-u[i]/0.00001;
        d[i]=-u[i];
    }
    else
    {
        if(v[i]>0.00001)
        d[i]=-u[i]/v[i];
        else
//        d[i]=-u[i]/0.00001;
        d[i]=-u[i];
    }
    b=-p*d;

    return 1;
}


int SysMatrixCanonical(matrix& a, matrix& p, vector& d)
{
    double eps=0.0000001;
    if(a.mrow()!=a.mcol())
    {
        cout<< Merrorstring[EMAT_NOTSQUARE];
        throw Merrorstring[EMAT_NOTSQUARE];
        return 0;
    }
try
{
    double dmin,d2;
    int i,j,k,m;
    m=a.mrow();
    p=eye(a.mrow(),a.mrow());
    d.initial(a.mrow());
    for(i=0;i<m-1;i++)
    {
        dmin=fabs(a(i,i));
        j=i;
        for(k=i+1;k<m;k++)
        {
            d2=fabs(a(i,k));
            if(dmin<d2)
            {
                dmin=d2;
                j=k;
            }
        }
        if(dmin<=eps)
        continue;
        if(j!=i)
        {
            d2=a(i,i)*a(i,j);
            if(d2>=0)
            {
                for(k=i;k<m;k++)
                {
                        a(i,k)+=a(j,k);
                }
                for(k=i;k<m;k++)
                {
                        a(k,i)+=a(k,j);
                }
                for(k=0;k<m;k++)
                {
                        p(k,i)+=p(k,j);
                }
            }
            else
            {
                for(k=i;k<m;k++)
                {
                        a(i,k)-=a(j,k);
                }
                for(k=i;k<m;k++)
                {
                        a(k,i)-=a(k,j);
                }
                for(k=0;k<m;k++)
                {
                        p(k,i)-=p(k,j);
                }
            }
        }
        for(k=i+1;k<m;k++)
        {
            d2=a(i,k);
            for(j=i;j<m;j++)
            a(j,k)=a(j,k)-d2*a(j,i)/a(i,i);
            for(j=0;j<m;j++)
            p(j,k)=p(j,k)-d2*p(j,i)/a(i,i);
        }
        for(k=i+1;k<m;k++)
        {
            d2=a(k,i);
            for(j=i;j<m;j++)
            a(k,j)=a(k,j)-d2*a(i,j)/a(i,i);
        }
    }
    for(i=0;i<a.mrow();i++)
    {
        d2=p(i,0);
        dmin=d2*d2;
        for(j=1;j<a.mrow();j++)
        {
                d2=p(i,j);
                dmin+=(d2*d2);
        }
        a(i,i)*=dmin;
        dmin=sqrt(dmin);
        p[i]/=dmin;
        d[i]=a(i,i);
    }
}
catch(char* str)
{
    cout<<"SysEig="<<str<<endl;
}
catch(...)
{
    cout<<"SysEig other exception"<<endl;
}
    return 1;
}

int CholeskyDecompose(matrix& A, matrix& L, vector& D,vector& E)
{

    int mr=A.mrow();
try
{
    L=eye(mr,mr);
    D.initial(mr);
    E.initial(mr);
    double d1=fabs(A(0,0));
    for(int k=1;k<mr;k++)
    {
        double data=fabs(A(k,k));
        if(d1<data)
        d1=data;
    }
    double d2=0;
    for(k=0;k<mr-1;k++)
    {
        for(int i=k+1;i<mr;i++)
        {
            double data=fabs(A(k,i));
            if(d2<data)
            d2=data;
        }
    }
    d2/=mr;
    if(d1<d2)
    d1=d2;
    d2=0.0000001;
    double dj=fabs(A(0,0));
    if(dj<d2)
    dj=d2;
    double dmax=0;
    for(k=1;k<mr;k++)
    {
        double data=fabs(A(k,0));
        if(dmax<data)
        dmax=data;
    }
    dmax=dmax*dmax/d1;
    if(dj<dmax)
    dj=dmax;
    for(k=1;k<mr;k++)
    L(k,0)=A(k,0)/dj;
    E[0]=dj-A(0,0);
	D[0]=dj;
    for(int j=1;j<mr;j++)
    {
        dmax=0;
		double data;
        for(int i=j+1;i<mr;i++)
        {
            data=0;
            for(k=0;k<j;k++)
            data+=A(j,k)*L(i,k);
            A(i,j)-=data;
			data=A(i,j);
            data=fabs(A(i,j));
            if(dmax<data)
            dmax=data;
        }
        dmax=dmax*dmax/d1;
		data=0;
		for(k=0;k<j;k++)
			data+=A(j,k)*L(j,k);
		dj=fabs(A(j,j)-data);
        if(dj<d2)
        dj=d2;
        if(dj<dmax)
        dj=dmax;
        for(k=j+1;k<mr;k++)
        L(k,j)=A(k,j)/dj;
        E[j]=dj+data-A(j,j);
		D[j]=dj;
    }
}
catch(char* str)
{
    cout<<"SysEig="<<str<<endl;
}
catch(...)
{
    cout<<"SysEig other exception"<<endl;
}
	return 1;
}

int SolveNewtonDir3(matrix& a,vector& b)
{
    vector v;
//    matrix L;
    v=copy(b);
//    double delta,beta,d1,d2;
    int mr,mc;
    mr=a.mrow();
    mc=a.mcol();
    if(mr!=mc)
    return 0;
    if(mr!=b())
    return 0;
    matrix p;
    vector e;
    vector d;
    CholeskyDecompose(a,p,v,e);
    d.initial(mr);
    for(int i=0;i<mr;i++)
    {
        double data=-b[i];
        for(int k=0;k<i;k++)
        data-=d[i]*p(i,k);
        d[i]=data;
    }
    for(i=0;i<mr;i++)
    {
        d[i]/=v[i];
    }
    vector b1;
    b1=copy(b);
    for(i=mr-1;i>=0;i--)
    {
        double data=0;
        for(int k=mr-1;k>i;k--)
        data+=b1[i]*p(k,i);
        b1[i]=(d[i]-data);
    }
    int k=-1;
    double d2=0;
    for(i=0;i<mr;i++)
    {
        double data=v[i]-e[i];
        if(data<d2)
        {
            d2=data;
            k=i;
        }
    }
    if(k>=0)
    {
        d=0;
        d[k]=1;
        for(i=k-1;i>=0;i--)
        {
            d[i]=0;
            for(int s=i+1;s<=k;s++)
            d[i]-=d[s]*p[s][i];
        }
        d2=dot(d,b);
        if(d2>0)
        {
            d=-d;
            d2=-d2;
        }
        d2+=v[k]-e[k];
        for(i=0;i<k;i++)
        d2-=d[i]*d[i]*e[i];
    }
    else
    d2=0;
    d2*=2;
    double d1=dot(b1,b);
    if(d1>0)
    {
        b1=-b1;
    }
    d1=fabs(d1);
    d2=fabs(d2);
    if(d1<=d2)
    {
        b=copy(d);
        cout<<1<<endl;;
    }
    else
    {
        b=copy(b1);
    }
    return 1;
}

void fmindir(double& g1,double& g2, double a,double b,double c)
{
    double delta=c*c-a*b;
    double d1,d2;
    if(delta<0)
    {
        d1=(c*g2-b*g1)/delta;
        d2=(c*g1-a*g2)/delta;
        g1=d1;g2=d2;
    }
    else
    {
        double data=a-b;
        double p11,p21,p12,p22;
        data=data*data+4*c*c;
        d1=(a+b+sqrt(data))/2;
        d2=(a+b-sqrt(data))/2;
        data=b-d1;
        p11=d1-b;
        p21=c;
        p12=-c;
        p22=p11;
        double h1,h2;
        data=p11*p11+p21*p21;
        h1=-(g1*p11+g2*p21)/data;
        h2=(g1*p12+g2*p22)/data;
        if(d1>0.00000001)
        h1/=d1;
        else
        h1/=0.00000001;
        if(d2<-0.00000001)
        h2/=d2;
        else
        h2/=-0.00000001;
        g1=h1*p11+h2*p12;
        g2=h1*p21+h2*p22;
    }
}

int SolveMNewtonDir(matrix& a,vector& b)
{
    matrix q;
    vector d,c;
try
{
    eastrq(a, q, d, c);
    double eps=DBL_EPSILON;
    int mr=d();
    int* I=new int[mr];
    vector dv(mr);
    double delta;
    for(int k=0;k<mr;k++)
    {
        double data=fabs(d[k]);
        double alfa;
        double beta;
        double gama;
        if( k<(mr-2) )
        {
            if(fabs(c[k])<=eps)
            {
                c[k]=0;
                I[k]=1;
            }
            else
            {
                delta=d[k]*d[k+1]-c[k]*c[k];
                if( data>=fabs(delta) )
                {
                    c[k]/=d[k];
                    I[k]=1;
                    d[k+1]-=c[k]*c[k];
                }
                else
                {
                    dv[k]=c[k];