NumericalTools.cpp

00001 #include "NumericalTools.h"
00002 
00003 extern StreamOutput cout;
00004 
00008 NumericalTools::NumericalTools()
00009 {
00010 }
00011 
00015 NumericalTools::~NumericalTools()
00016 {
00017 }
00018 
00033 void NumericalTools::RidgeRegressionMultisolutionSinglecallGELSS ( REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, bool isRowMajor )
00034 {
00035     int oneInt=1, rank, info;
00036     REAL precision = 1e-15;
00037 
00038     // the new one
00039     int nn = n + m;
00040     int lWork = -1;
00041     REAL workSize[1];
00042     REAL* AA = new REAL[nn*m];
00043     REAL* bb = new REAL[nn*k];
00044     REAL* eigenvals = new REAL[m];
00045 
00046     GELSS ( &nn, &m, &k, AA, &nn, bb, &nn, eigenvals, &precision, &rank, workSize, &lWork, &info );
00047 
00048     lWork = ( int ) ( workSize[0]+0.5 );
00049     REAL work[lWork];
00050 
00051     for ( int i=0;i<n;i++ )
00052     {
00053         for ( int j=0;j<m;j++ )
00054             AA[i + j*nn] = A[i + j*n];
00055         for ( int j=0;j<k;j++ )
00056             bb[i + j*nn] = b[i + j*n];
00057     }
00058 
00059     // add regularizer matrix to the bottom of A
00060     //
00061     // A = [1 2 3]  b = [1]
00062     //     [3 4 5]      [2]
00063     //     [.. ..]      [.]
00064     // lam*[1 0 0]      [0]  -> solving this extended equation system is
00065     // lam*[0 1 0]      [0]  -> equivalent to ridge regression !
00066     // lam*[0 0 1]      [0]  -> (Paper from Y.Koren 2007 used that solving schema)
00067     REAL lam = sqrt ( lambda* ( REAL ) n );
00068     for ( int i=0;i<m;i++ )
00069     {
00070         for ( int j=0;j<m;j++ )
00071         {
00072             if ( i == j )
00073                 AA[j + i*nn + n] = lam;
00074             else
00075                 AA[j + i*nn + n] = 0.0;
00076         }
00077         bb[i + n] = 0.0;
00078     }
00079 
00080     GELSS ( &nn, &m, &k, AA, &nn, bb, &nn, eigenvals, &precision, &rank, work, &lWork, &info );
00081 
00082     for ( int i=0;i<m;i++ )
00083         for ( int j=0;j<k;j++ )
00084             x[i+j*k] = bb[i+j*nn];
00085 
00086     if ( AA )
00087         delete[] AA;
00088     AA = 0;
00089     if ( bb )
00090         delete[] bb;
00091     bb = 0;
00092     if ( eigenvals )
00093         delete[] eigenvals;
00094     eigenvals = 0;
00095 }
00096 
00110 void NumericalTools::RidgeRegressionMultisolutionMulticall ( REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, int alloc, bool isRowMajor )
00111 {
00112     static REAL *AA;
00113     static REAL *Ab;
00114     static REAL *AbTrans;
00115     static REAL *eigenvals;
00116     static REAL *work;
00117     static int *iwork;
00118     int oneInt=1, rank, info;
00119     REAL precision = 1e-15;
00120     static int lWork=-1;
00121     static REAL workSize[1];
00122 
00123     if ( alloc == ALLOC_MEM )
00124     {
00125         AA = new REAL[m*m];
00126         Ab = new REAL[m*k];
00127         AbTrans = new REAL[m*k];
00128         eigenvals = new REAL[m];
00129         lWork = -1;
00130         GELSS ( &m, &m, &k, AA, &m, Ab, &m, eigenvals, &precision, &rank, workSize, &lWork, &info );
00131 
00132         lWork = ( int ) ( workSize[0]+0.5 );
00133         work = new REAL[lWork];
00134 
00135         return;
00136     }
00137 
00138     if ( alloc == FREE_MEM )
00139     {
00140         if ( work )
00141             delete[] work;
00142         work = 0;
00143         if ( AA )
00144             delete[] AA;
00145         AA = 0;
00146         if ( Ab )
00147             delete[] Ab;
00148         Ab = 0;
00149         if ( AbTrans )
00150             delete[] AbTrans;
00151         AbTrans = 0;
00152         if ( eigenvals )
00153             delete[] eigenvals;
00154         eigenvals = 0;
00155 
00156         return;
00157     }
00158 
00159     if ( isRowMajor )
00160         CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, AA, m ); // AA = A'A
00161     else
00162         CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, AA, m ); // AA = A'A
00163 
00164     for ( int i=0; i < m; i++ )
00165         AA[i+i*m] += ( REAL ) n*lambda;  // A'A + lambda*I
00166 
00167     if ( isRowMajor )
00168         CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, k , n, 1.0, A, m, b, k, 0.0, Ab, k ); // Ab = A'b
00169     else
00170         CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b
00171 
00172     if ( isRowMajor )
00173     {
00174         for ( int i=0;i<k;i++ )
00175             for ( int j=0;j<m;j++ )
00176                 AbTrans[j+i*m] = Ab[i+j*k];  // copy to colMajor
00177     }
00178     else
00179     {
00180         for ( int i=0;i<k;i++ )
00181             for ( int j=0;j<m;j++ )
00182                 AbTrans[i+j*k] = Ab[i+j*k];  // copy to colMajor
00183     }
00184 
00185     GELSS ( &m, &m, &k, AA, &m, AbTrans, &m, eigenvals, &precision, &rank, work, &lWork, &info );
00186 
00187     if ( info != 0 )
00188         cout<<"error equation solver failed to converge: "<<info<<endl;
00189 
00190     if ( isRowMajor )
00191     {
00192         for ( int i=0;i<m;i++ ) // all features
00193             for ( int j=0;j<k;j++ ) // all classes (#outupts)
00194                 x[j+i*k] = AbTrans[i+j*m];
00195     }
00196     else
00197     {
00198         for ( int i=0;i<m;i++ ) // all features
00199             for ( int j=0;j<k;j++ ) // all classes (#outupts)
00200                 x[i+j*m] = AbTrans[i+j*m];
00201     }
00202 
00203 }
00204 
00218 void NumericalTools::RidgeRegressionMultisolutionSinglecall ( REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, bool isRowMajor, double* rideModifier )
00219 {
00220     REAL *AA;
00221     REAL *Ab;
00222     REAL *AbTrans;
00223     REAL *eigenvals;
00224     //REAL *work;
00225     int *iwork;
00226     int oneInt=1, rank, info;
00227     REAL precision = 1e-15;
00228     int lWork=-1;
00229     REAL workSize[1];
00230 
00231     AA = new REAL[m*m];
00232     Ab = new REAL[m*k];
00233     AbTrans = new REAL[m*k];
00234     eigenvals = new REAL[m];
00235     lWork = -1;
00236 
00237     // get size of work array
00238     //GELSS(&m, &m, &k, AA, &m, Ab, &m, eigenvals, &precision, &rank, workSize, &lWork, &info);
00239     //lWork = (int)(workSize[0]+0.5);
00240     //work = new REAL[lWork];
00241 
00242     if ( isRowMajor )
00243         CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, AA, m ); // AA = A'A
00244     else
00245         CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, AA, m ); // AA = A'A
00246 
00247     if ( rideModifier )
00248         for ( int i=0; i < m; i++ )
00249             AA[i+i*m] += ( REAL ) n*lambda*rideModifier[i];  // A'A + lambda*I*modifier
00250     else
00251         for ( int i=0; i < m; i++ )
00252             AA[i+i*m] += ( REAL ) n*lambda;  // A'A + lambda*I
00253 
00254     if ( isRowMajor )
00255         CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, k , n, 1.0, A, m, b, k, 0.0, Ab, k ); // Ab = A'b
00256     else
00257         CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b
00258 
00259     if ( isRowMajor )
00260     {
00261         for ( int i=0;i<k;i++ )
00262             for ( int j=0;j<m;j++ )
00263                 AbTrans[j+i*m] = Ab[i+j*k];  // copy to colMajor
00264     }
00265     else
00266     {
00267         for ( int i=0;i<k;i++ )
00268             for ( int j=0;j<m;j++ )
00269                 AbTrans[i+j*k] = Ab[i+j*k];  // copy to colMajor
00270     }
00271 
00272     int lda = m, nrhs = k, ldb = m;
00273 
00274     //Cholesky factorization of a symmetric (Hermitian) positive-definite matrix
00275     LAPACK_POTRF ( "U", &m, AA, &lda, &info );
00276 
00277     // Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite matrix
00278     LAPACK_POTRS ( "U", &m, &nrhs, AA, &lda, AbTrans, &ldb, &info );  // Y=solution vector
00279 
00280 
00281     // slower solver
00282     // Computes the minimum-norm solution to a linear least squares problem using the singular value decomposition of A
00283     //REAL *work;
00284     // get size of work array
00285     //GELSS(&m, &m, &k, AA, &m, Ab, &m, eigenvals, &precision, &rank, workSize, &lWork, &info);
00286     //lWork = (int)(workSize[0]+0.5);
00287     //work = new REAL[lWork];
00288     //GELSS(&m, &m, &k, AA, &m, AbTrans, &m, eigenvals, &precision, &rank, work, &lWork, &info);
00289     //delete[] work;
00290 
00291     if ( info != 0 )
00292         cout<<"error equation solver failed to converge: "<<info<<endl;
00293 
00294     if ( isRowMajor )
00295     {
00296         for ( int i=0;i<m;i++ ) // all features
00297             for ( int j=0;j<k;j++ ) // all classes (#outupts)
00298                 x[j+i*k] = AbTrans[i+j*m];
00299     }
00300     else
00301     {
00302         for ( int i=0;i<m;i++ ) // all features
00303             for ( int j=0;j<k;j++ ) // all classes (#outupts)
00304                 x[i+j*m] = AbTrans[i+j*m];
00305     }
00306 
00307     if ( AA )
00308         delete[] AA;
00309     AA = 0;
00310     if ( Ab )
00311         delete[] Ab;
00312     Ab = 0;
00313     if ( AbTrans )
00314         delete[] AbTrans;
00315     AbTrans = 0;
00316     if ( eigenvals )
00317         delete[] eigenvals;
00318     eigenvals = 0;
00319 }
00320 
00337 void NumericalTools::LogisticRegressionMultisolutionSinglecall ( REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, REAL offset, REAL scale, bool isRowMajor )
00338 {
00339     /* logistic regression - MATLAB
00340     function x = logistic4(A, y, K, ridge, epsilon)
00341     N = size(A,1);
00342     M = size(A,2);
00343 
00344     pi = ones(N,K) * 1/K;
00345     eta = log(pi);
00346 
00347     kron1 = repmat(1:K-1,M,1);
00348     kron2 = repmat((1:M)',1,K-1);
00349 
00350     x = zeros(M*(K-1),1);
00351     xold = x;
00352 
00353     for epoch=0:100
00354         mu = pi;
00355         XWX = zeros(M*(K-1),M*(K-1));
00356         XWZ = zeros(M*(K-1),1);
00357         for i=1:N  % over all training samples
00358             W = diag(mu(i,:)) - mu(i,:)'*pi(i,:);
00359             Z = eta(i,:)*W + (y(i,:) - mu(i,:));
00360             astar = A(i,:);
00361             XWX = XWX + W(kron1,kron1) .* (astar(1,kron2)'*astar(1,kron2));
00362             XWZ = XWZ + Z(1,kron1)' .* astar(1,kron2)';
00363         end
00364         x = (XWX + eye(M*(K-1))*rigde*M*(K-1)) \ XWZ;
00365         x = reshape(x,M,K-1);
00366         disp(['x:' num2str(x(:)')]);
00367         if sum((x(:)-xold(:)).^2)/(N*(K-1)) < epsilon
00368             return;
00369         end
00370         xold = x;
00371         etaOld = eta;
00372         eta = A*x;
00373         eta = [eta zeros(N,1)];
00374         
00375         for backstep = 0:10
00376             pi = exp(eta);
00377             pi = pi ./ repmat(sum(pi,2),1,K);
00378             if all(pi(:) > 1.8e-12)
00379                 break;
00380             elseif backstep < 10
00381                 eta = etaOld + (eta - etaOld)/5;
00382             else
00383                 pi = max(pi,1.8e-12);
00384                 pi = pi ./ repmat(sum(pi,2),1,K);
00385                 eta = log(pi);
00386                 break;
00387             end
00388         end
00389     end
00390     */
00391     
00392     REAL tolerance = 1e-5; //sizeof(REAL)==4 ? 6.4155e-06 : 1.8190e-12;
00393     
00394     REAL* xOld = new REAL[m*(k-1)];
00395     REAL* pi = new REAL[n*k];
00396     REAL* eta = new REAL[n*k];
00397     REAL* etaTmp = new REAL[n*(k-1)];
00398     REAL* etaOld = new REAL[n*k];
00399     REAL* mu = new REAL[n*k];
00400     REAL* XWX = new REAL[m*(k-1) * m*(k-1)];
00401     REAL* XWZ = new REAL[m*(k-1)];
00402     REAL* W = new REAL[k*k];
00403     REAL* Z = new REAL[k];
00404     for(int i=0;i<n*k;i++)
00405     {
00406         pi[i] = 1.0 / (REAL)k;
00407         eta[i] = log(pi[i]);
00408     }
00409     for(int i=0;i<m*(k-1);i++)
00410         xOld[i] = 0.0;
00411     
00412     for(int epoch=0;epoch<100;epoch++)  // IRLS - iterative reweighted least squares loop
00413     {
00414         for(int i=0;i<n*k;i++)
00415             mu[i] = pi[i];
00416         for(int i=0;i<m*(k-1) * m*(k-1);i++)
00417             XWX[i] = 0.0;
00418         for(int i=0;i<m*(k-1);i++)
00419             XWZ[i] = 0.0;
00420         
00421         for(int e=0;e<n;e++) // loop over training examples
00422         {
00423             REAL *muPtr = mu + e*k, *piPtr = pi + e*k, *etaPtr = eta + e*k, *bPtr = b + e*k, *Aptr = A + e*m;
00424             for(int i=0;i<k;i++)
00425                 for(int j=0;j<k;j++)
00426                     W[i*k+j] = (i==j? muPtr[i] : 0.0) - muPtr[i] * piPtr[j];
00427             for(int i=0;i<k;i++)
00428             {
00429                 REAL sum = 0.0;
00430                 for(int j=0;j<k;j++)
00431                     sum += etaPtr[j] * W[i+j*k];
00432                 Z[i] = sum + (bPtr[i] - muPtr[i]);
00433             }
00434             // this is too slow code
00435             for(int i=0;i<m*(k-1);i++)
00436                 for(int j=0;j<m*(k-1);j++)
00437                     XWX[i*m*(k-1) + j] += W[(i/m)*k+j/m] * Aptr[i%m] * Aptr[j%m];
00438             for(int i=0;i<k-1;i++)
00439                 for(int j=0;j<m;j++)
00440                     XWZ[i*m+j] += Z[i] * Aptr[j];
00441         }
00442         
00443         for(int i=0;i<m*(k-1);i++)
00444             XWX[i+i*m*(k-1)] += lambda*(REAL)m*(k-1);
00445         
00446         int lda = m*(k-1), nrhs = 1, ldb = m*(k-1), info;
00447         //Cholesky factorization of a symmetric (Hermitian) positive-definite matrix
00448         LAPACK_POTRF ( "U", &lda, XWX, &lda, &info );
00449         // Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite matrix
00450         LAPACK_POTRS ( "U", &lda, &nrhs, XWX, &lda, XWZ, &ldb, &info );  // Ab=solution vector
00451         if ( info != 0 )
00452             cout<<"error equation solver failed to converge: "<<info<<endl;
00453         for(int i=0;i<k-1;i++)
00454             for(int j=0;j<m;j++)
00455                 x[i+j*(k-1)] = XWZ[j+i*m];
00456         
00457         //for(int i=0;i<m*(k-1);i++)
00458         //    cout<<x[i]<<" ";
00459         
00460         double sum0 = 0.0, sum1 = 0.0, normalizedDiff;
00461         for(int i=0;i<m*(k-1);i++)
00462         {
00463             sum0 += x[i] * x[i];
00464             sum1 += xOld[i] * xOld[i];
00465         }
00466         normalizedDiff = (sum0-sum1)/((double)m*(k-1));
00467         //cout<<" | "<<normalizedDiff<<endl;
00468         if(normalizedDiff < tolerance)  // break criteria
00469             return;
00470         
00471         for(int i=0;i<m*(k-1);i++)
00472             xOld[i] = x[i];
00473         
00474         for(int i=0;i<n*k;i++)
00475             etaOld[i] = eta[i];
00476         
00477         CBLAS_GEMM ( CblasRowMajor, CblasNoTrans, CblasNoTrans, n, k-1 , m, 1.0, A, m, x, k-1, 0.0, etaTmp, k-1 ); // etaTmp = A * x
00478         
00479         for(int i=0;i<n;i++)
00480         {
00481             for(int j=0;j<k-1;j++)
00482                 eta[j+i*k] = etaTmp[j+i*(k-1)];
00483             eta[k-1+i*k] = 0.0;
00484         }
00485         
00486         for(int backstep=0;backstep<11;backstep++)
00487         {
00488             for(int i=0;i<n*k;i++)
00489                 pi[i] = exp(eta[i]);
00490             
00491             for(int i=0;i<n;i++)
00492             {
00493                 REAL sum = 0.0;
00494                 REAL* ptr = pi + i*k;
00495                 for(int j=0;j<k;j++)
00496                     sum += ptr[j];
00497                 for(int j=0;j<k;j++)
00498                     ptr[j] /= sum;
00499             }
00500             
00501             uint sum = 0;
00502             for(int i=0;i<n*k;i++)
00503                 if(pi[i] > tolerance)
00504                     sum++;
00505             if(sum == n*k)
00506                 break;
00507             else if(backstep < 10)
00508             {
00509                 for(int i=0;i<n*k;i++)
00510                     eta[i] = etaOld[i] + (eta[i] - etaOld[i])*0.2;
00511             }
00512             else
00513             {
00514                 for(int i=0;i<n*k;i++)
00515                     if(pi[i] <= tolerance)
00516                         pi[i] = tolerance;
00517                 for(int i=0;i<n;i++)
00518                 {
00519                     REAL sum = 0.0;
00520                     REAL* ptr = pi + i*k;
00521                     for(int j=0;j<k;j++)
00522                         sum += ptr[j];
00523                     for(int j=0;j<k;j++)
00524                         pi[j] /= sum;
00525                 }
00526                 for(int i=0;i<n*k;i++)
00527                     eta[i] = log(pi[i]);
00528                 break;
00529             }
00530         }
00531     }
00532     
00533     delete[] xOld;
00534     delete[] pi;
00535     delete[] eta;
00536     delete[] etaTmp;
00537     delete[] etaOld;
00538     delete[] mu;
00539     delete[] XWX;
00540     delete[] XWZ;
00541     delete[] W;
00542     delete[] Z;
00543 }
00544 
00562 int NumericalTools::RidgeRegressionNonNegSinglecall ( REAL *A, REAL *b, REAL *x, int n, int m, REAL lamda, REAL eps, int maxIter, bool isRowMajor, int debug )
00563 {
00564     static REAL *r, *tmp, *ATA, *Ab;
00565 
00566     r = new REAL[m];
00567     tmp = new REAL[m];
00568     ATA = new REAL[m*m];
00569     Ab = new REAL[m];
00570 
00571     int nIter = 0;
00572     unsigned int i;
00573     REAL alpha = 0.0, tmp0, tmp1;
00574     REAL norm = 1.0;
00575 
00576     if ( lamda > 0.0 )
00577     {
00578         if ( isRowMajor )
00579             CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, ATA, m ); // AA = A'A
00580         else
00581             CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, ATA, m ); // AA = A'A
00582 
00583         if ( isRowMajor )
00584             CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, m, b, 1, 0.0, Ab, 1 ); // calc Ab = A'*B
00585         else
00586             CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b
00587 
00588         // add constant diagonal elements for regularization
00589         for ( i=0;i<m;i++ )
00590             ATA[i+i*m] += lamda * ( REAL ) n;
00591 
00592     }
00593     else
00594     {
00595         memcpy ( ( char* ) ATA, ( char* ) A, sizeof ( REAL ) *m*m );
00596         memcpy ( ( char* ) Ab, ( char* ) b, sizeof ( REAL ) *m );
00597     }
00598 
00599     // init solution
00600     for ( i=0;i<m;i++ )
00601         x[i] = r[i] = 1.0;
00602 
00603     while ( norm>eps )
00604     {
00605         CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, x, 1, 0.0, tmp, 1 );  // tmp = ATA*x
00606         for ( i=0;i<m;i++ )
00607             r[i] = Ab[i] - tmp[i];
00608         for ( i=0;i<m;i++ )
00609             if ( x[i]==0.0 &&  r[i]<0.0 )
00610                 r[i] = 0.0;
00611         CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, r, 1, 0.0, tmp, 1 );  // tmp = A*r
00612         tmp0 = tmp1 = 0.0;
00613         for ( i=0;i<m;i++ )
00614         {
00615             tmp0 += r[i]*r[i];   // nominator
00616             tmp1 += r[i]*tmp[i]; // denominator
00617         }
00618         if ( tmp0!=0.0 && tmp1!=0.0 )
00619             alpha = tmp0 / tmp1;  // alpha = r'r / r'Ar   (tmp=Ar)
00620         for ( i=0;i<m;i++ )
00621         {
00622             if ( r[i]<0.0 )
00623             {
00624                 tmp0 = -x[i]/r[i];
00625                 alpha = alpha<tmp0? alpha : tmp0;
00626             }
00627         }
00628         for ( i=0;i<m;i++ )
00629             x[i] += alpha * r[i];
00630         norm = 0.0;
00631         for ( i=0;i<m;i++ )
00632         {
00633             if ( x[i]<1e-9 )
00634                 x[i] = 0.0;  // x(x<1e-6) = 0
00635             norm += r[i]*r[i];
00636         }
00637         norm = sqrt ( norm );
00638         nIter++;
00639 
00640         // break criteria, if no convergence
00641         if ( nIter>=maxIter )
00642             break;
00643 
00644         if ( debug )
00645         {
00646             cout<<"norm="<<norm<<"  alpha="<<alpha<<endl;
00647             cout<<"r=";
00648             for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
00649             cout<<endl;
00650             cout<<"x=";
00651             for ( i=0;i<m;i++ ) cout<<x[i]<<" ";
00652             cout<<endl;
00653             cout<<"tmp=";
00654             for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
00655             cout<<endl;
00656             cout<<endl;
00657         }
00658     }
00659 
00660     // check the solution for nan of inf values
00661     for ( i=0;i<m;i++ )
00662         if ( isnan ( x[i] ) || isinf ( x[i] ) )
00663             cout<<"DIVERGENCE of nonNegSolver (nIter="<<nIter<<")"<<endl;
00664 
00665     if ( r )
00666         delete[] r;
00667     r = 0;
00668     if ( tmp )
00669         delete[] tmp;
00670     tmp = 0;
00671     if ( ATA )
00672         delete[] ATA;
00673     ATA = 0;
00674     if ( Ab )
00675         delete[] Ab;
00676     Ab = 0;
00677 
00678     return nIter;
00679 }
00680 
00700 int NumericalTools::RidgeRegressionNonNegmulticall ( REAL *A, REAL *b, REAL *x, int n, int m, REAL lamda, REAL eps, int maxIter, int alloc, bool isRowMajor, int debug )
00701 {
00702     static REAL *r, *tmp, *ATA, *Ab;
00703 
00704     if ( alloc == ALLOC_MEM )
00705     {
00706         cout<<"Allocate memory in nonNeg solver"<<endl;
00707         r = new REAL[m];
00708         tmp = new REAL[m];
00709         ATA = new REAL[m*m];
00710         Ab = new REAL[m];
00711         return -1;
00712     }
00713     if ( alloc == FREE_MEM )
00714     {
00715         cout<<"Free memory in nonNeg solver"<<endl;
00716         if ( r )
00717             delete[] r;
00718         r = 0;
00719         if ( tmp )
00720             delete[] tmp;
00721         tmp = 0;
00722         if ( ATA )
00723             delete[] ATA;
00724         ATA = 0;
00725         if ( Ab )
00726             delete[] Ab;
00727         Ab = 0;
00728         return -1;
00729     }
00730 
00731     int nIter = 0;
00732     unsigned int i;
00733     REAL alpha = 0.0, tmp0, tmp1;
00734     REAL norm = 1.0;
00735 
00736     if ( lamda > 0.0 )
00737     {
00738         if ( isRowMajor )
00739             CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, ATA, m ); // AA = A'A
00740         else
00741             CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, ATA, m ); // AA = A'A
00742 
00743         if ( isRowMajor )
00744             CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, m, b, 1, 0.0, Ab, 1 ); // calc Ab = A'*B
00745         else
00746             CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b
00747 
00748         // add constant diagonal elements for regularization
00749         for ( i=0;i<m;i++ )
00750             ATA[i+i*m] += lamda * ( REAL ) n;
00751 
00752     }
00753     else
00754     {
00755         memcpy ( ( char* ) ATA, ( char* ) A, sizeof ( REAL ) *m*m );
00756         memcpy ( ( char* ) Ab, ( char* ) b, sizeof ( REAL ) *m );
00757     }
00758 
00759     // init solution
00760     for ( i=0;i<m;i++ )
00761         x[i] = r[i] = 1.0;
00762 
00763     while ( norm>eps )
00764     {
00765         CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, x, 1, 0.0, tmp, 1 );  // tmp = ATA*x
00766         for ( i=0;i<m;i++ )
00767             r[i] = Ab[i] - tmp[i];
00768         for ( i=0;i<m;i++ )
00769             if ( x[i]==0.0 &&  r[i]<0.0 )
00770                 r[i] = 0.0;
00771         CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, r, 1, 0.0, tmp, 1 );  // tmp = A*r
00772         tmp0 = tmp1 = 0.0;
00773         for ( i=0;i<m;i++ )
00774         {
00775             tmp0 += r[i]*r[i];   // nominator
00776             tmp1 += r[i]*tmp[i]; // denominator
00777         }
00778         if ( tmp0!=0.0 && tmp1!=0.0 )
00779             alpha = tmp0 / tmp1;  // alpha = r'r / r'Ar   (tmp=Ar)
00780         for ( i=0;i<m;i++ )
00781         {
00782             if ( r[i]<0.0 )
00783             {
00784                 tmp0 = -x[i]/r[i];
00785                 alpha = alpha<tmp0? alpha : tmp0;
00786             }
00787         }
00788         for ( i=0;i<m;i++ )
00789             x[i] += alpha * r[i];
00790         norm = 0.0;
00791         for ( i=0;i<m;i++ )
00792         {
00793             if ( x[i]<1e-9 )
00794                 x[i] = 0.0;  // x(x<1e-6) = 0
00795             norm += r[i]*r[i];
00796         }
00797         norm = sqrt ( norm );
00798         nIter++;
00799 
00800         // break criteria, if no convergence
00801         if ( nIter>=maxIter )
00802             break;
00803 
00804         if ( debug )
00805         {
00806             cout<<"norm="<<norm<<"  alpha="<<alpha<<endl;
00807             cout<<"r=";
00808             for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
00809             cout<<endl;
00810             cout<<"x=";
00811             for ( i=0;i<m;i++ ) cout<<x[i]<<" ";
00812             cout<<endl;
00813             cout<<"tmp=";
00814             for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
00815             cout<<endl;
00816             cout<<endl;
00817         }
00818     }
00819 
00820     // check the solution for nan of inf values
00821     for ( i=0;i<m;i++ )
00822         if ( isnan ( x[i] ) || isinf ( x[i] ) )
00823             cout<<"DIVERGENCE of nonNegSolver (nIter="<<nIter<<")"<<endl;
00824 
00825     return nIter;
00826 }
00827 
00835 REAL NumericalTools::getNormRandomNumber ( REAL mean, REAL std )
00836 {
00837     REAL x1, x2, w, y1;
00838 
00839     do
00840     {
00841         x1 = 2.0 * ( ( REAL ) rand() / ( REAL ) RAND_MAX ) - 1.0;
00842         x2 = 2.0 * ( ( REAL ) rand() / ( REAL ) RAND_MAX ) - 1.0;
00843         w = x1 * x1 + x2 * x2;
00844     }
00845     while ( w >= 1.0 );
00846 
00847     w = sqrt ( ( -2.0 * log ( w ) ) / w );
00848     y1 = x1 * w;
00849 
00850     y1 *= std;
00851     y1 += mean;
00852 
00853     return y1;
00854 }
00855 
00863 REAL NumericalTools::getUniformRandomNumber ( REAL min, REAL max )
00864 {
00865     REAL  x = ( REAL ) rand() / ( REAL ) RAND_MAX;
00866 
00867     x *= ( max-min );
00868     x = min + x;
00869 
00870     return x;
00871 }
00872 
00880 REAL NumericalTools::clipValue ( REAL value, REAL min, REAL max )
00881 {
00882     if ( value < min )
00883         value = min;
00884     if ( value > max )
00885         value = max;
00886 
00887     return value;
00888 }
Generated on Tue Jan 26 09:20:59 2010 for ELF by  doxygen 1.5.8