NumericalTools Class Reference

#include <NumericalTools.h>

Inheritance diagram for NumericalTools:

List of all members.

Public Member Functions
	NumericalTools ()
	~NumericalTools ()
int	RidgeRegressionNonNegSinglecall (REAL *A, REAL *b, REAL *x, int n, int m, REAL lamda, REAL eps, int maxIter, bool isRowMajor=true, int debug=0)
Static Public Member Functions
static void	RidgeRegressionMultisolutionSinglecallGELSS (REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, bool isRowMajor=false)
static void	RidgeRegressionMultisolutionSinglecall (REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, bool isRowMajor, double *rideModifier=0)
static void	RidgeRegressionMultisolutionMulticall (REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, int alloc=0, bool isRowMajor=false)
static int	RidgeRegressionNonNegmulticall (REAL *A, REAL *b, REAL *x, int n, int m, REAL lamda, REAL eps, int maxIter, int alloc=0, bool isRowMajor=true, int debug=0)
static void	LogisticRegressionMultisolutionSinglecall (REAL *A, REAL *b, REAL *x, int n, int m, int k, REAL lambda, REAL offset=0.0, REAL scale=1.0, bool isRowMajor=true)
static REAL	getNormRandomNumber (REAL mean, REAL std)
static REAL	getUniformRandomNumber (REAL min, REAL max)
static REAL	clipValue (REAL value, REAL min, REAL max)

---

Detailed Description

A collection of numerical tools that are used in this framework

• linear equation solvers (for multiple targets, non-negative solver)
• logistic equation solver (experimental)
• helper fuctions

Definition at line 30 of file NumericalTools.h.

---

Constructor & Destructor Documentation

NumericalTools::NumericalTools

(

)

Constructor

Definition at line 8 of file NumericalTools.cpp.

{
}

NumericalTools::~NumericalTools

(

)

Destructor

Definition at line 15 of file NumericalTools.cpp.

{
}

---

Member Function Documentation

REAL NumericalTools::clipValue	(	REAL	value,
		REAL	min,
		REAL	max
	)			[static]

Clip a value with lower and upper bound

Parameters:

	min	Min value
	max	Max value

Returns:

Clipped value

Definition at line 880 of file NumericalTools.cpp.

{
    if ( value < min )
        value = min;
    if ( value > max )
        value = max;

    return value;
}

REAL NumericalTools::getNormRandomNumber	(	REAL	mean,
		REAL	std
	)			[static]

Returns a gaussian random number

Parameters:

	mean	Mean value
	std	Standard deviation

Returns:

The random value

Definition at line 835 of file NumericalTools.cpp.

{
    REAL x1, x2, w, y1;

    do
    {
        x1 = 2.0 * ( ( REAL ) rand() / ( REAL ) RAND_MAX ) - 1.0;
        x2 = 2.0 * ( ( REAL ) rand() / ( REAL ) RAND_MAX ) - 1.0;
        w = x1 * x1 + x2 * x2;
    }
    while ( w >= 1.0 );

    w = sqrt ( ( -2.0 * log ( w ) ) / w );
    y1 = x1 * w;

    y1 *= std;
    y1 += mean;

    return y1;
}

REAL NumericalTools::getUniformRandomNumber	(	REAL	min,
		REAL	max
	)			[static]

Returns a uniform random number

Parameters:

	min	Min value
	max	Max value

Returns:

Random value

Definition at line 863 of file NumericalTools.cpp.

{
    REAL  x = ( REAL ) rand() / ( REAL ) RAND_MAX;

    x *= ( max-min );
    x = min + x;

    return x;
}

void NumericalTools::LogisticRegressionMultisolutionSinglecall	(	REAL *	A,
		REAL *	b,
		REAL *	x,
		int	n,
		int	m,
		int	k,
		REAL	lambda,
		REAL	offset = 0.0,
		REAL	scale = 1.0,
		bool	isRowMajor = true
	)			[static]

Logistic Regression with multiple target vectors

Use this function, if you have to solve multiple logistic equations. The target vector b gets transformed by: bInternal = scale*b+offset;

A ... n x m matrix b ... n x k target vector (0s or 1s) x ... m x k solution vector

This is a 1:1 copy from the matlab function "mnrfit" -> Fit a nominal or ordinal multinomial regression model The mnrfit can be found under $MATLAB/toolbox/stats/mnrfit.m

Definition at line 337 of file NumericalTools.cpp.

{
    /* logistic regression - MATLAB
    function x = logistic4(A, y, K, ridge, epsilon)
    N = size(A,1);
    M = size(A,2);

    pi = ones(N,K) * 1/K;
    eta = log(pi);

    kron1 = repmat(1:K-1,M,1);
    kron2 = repmat((1:M)',1,K-1);

    x = zeros(M*(K-1),1);
    xold = x;

    for epoch=0:100
        mu = pi;
        XWX = zeros(M*(K-1),M*(K-1));
        XWZ = zeros(M*(K-1),1);
        for i=1:N  % over all training samples
            W = diag(mu(i,:)) - mu(i,:)'*pi(i,:);
            Z = eta(i,:)*W + (y(i,:) - mu(i,:));
            astar = A(i,:);
            XWX = XWX + W(kron1,kron1) .* (astar(1,kron2)'*astar(1,kron2));
            XWZ = XWZ + Z(1,kron1)' .* astar(1,kron2)';
        end
        x = (XWX + eye(M*(K-1))*rigde*M*(K-1)) \ XWZ;
        x = reshape(x,M,K-1);
        disp(['x:' num2str(x(:)')]);
        if sum((x(:)-xold(:)).^2)/(N*(K-1)) < epsilon
            return;
        end
        xold = x;
        etaOld = eta;
        eta = A*x;
        eta = [eta zeros(N,1)];
        
        for backstep = 0:10
            pi = exp(eta);
            pi = pi ./ repmat(sum(pi,2),1,K);
            if all(pi(:) > 1.8e-12)
                break;
            elseif backstep < 10
                eta = etaOld + (eta - etaOld)/5;
            else
                pi = max(pi,1.8e-12);
                pi = pi ./ repmat(sum(pi,2),1,K);
                eta = log(pi);
                break;
            end
        end
    end
    */
    
    REAL tolerance = 1e-5; //sizeof(REAL)==4 ? 6.4155e-06 : 1.8190e-12;
    
    REAL* xOld = new REAL[m*(k-1)];
    REAL* pi = new REAL[n*k];
    REAL* eta = new REAL[n*k];
    REAL* etaTmp = new REAL[n*(k-1)];
    REAL* etaOld = new REAL[n*k];
    REAL* mu = new REAL[n*k];
    REAL* XWX = new REAL[m*(k-1) * m*(k-1)];
    REAL* XWZ = new REAL[m*(k-1)];
    REAL* W = new REAL[k*k];
    REAL* Z = new REAL[k];
    for(int i=0;i<n*k;i++)
    {
        pi[i] = 1.0 / (REAL)k;
        eta[i] = log(pi[i]);
    }
    for(int i=0;i<m*(k-1);i++)
        xOld[i] = 0.0;
    
    for(int epoch=0;epoch<100;epoch++)  // IRLS - iterative reweighted least squares loop
    {
        for(int i=0;i<n*k;i++)
            mu[i] = pi[i];
        for(int i=0;i<m*(k-1) * m*(k-1);i++)
            XWX[i] = 0.0;
        for(int i=0;i<m*(k-1);i++)
            XWZ[i] = 0.0;
        
        for(int e=0;e<n;e++) // loop over training examples
        {
            REAL *muPtr = mu + e*k, *piPtr = pi + e*k, *etaPtr = eta + e*k, *bPtr = b + e*k, *Aptr = A + e*m;
            for(int i=0;i<k;i++)
                for(int j=0;j<k;j++)
                    W[i*k+j] = (i==j? muPtr[i] : 0.0) - muPtr[i] * piPtr[j];
            for(int i=0;i<k;i++)
            {
                REAL sum = 0.0;
                for(int j=0;j<k;j++)
                    sum += etaPtr[j] * W[i+j*k];
                Z[i] = sum + (bPtr[i] - muPtr[i]);
            }
            // this is too slow code
            for(int i=0;i<m*(k-1);i++)
                for(int j=0;j<m*(k-1);j++)
                    XWX[i*m*(k-1) + j] += W[(i/m)*k+j/m] * Aptr[i%m] * Aptr[j%m];
            for(int i=0;i<k-1;i++)
                for(int j=0;j<m;j++)
                    XWZ[i*m+j] += Z[i] * Aptr[j];
        }
        
        for(int i=0;i<m*(k-1);i++)
            XWX[i+i*m*(k-1)] += lambda*(REAL)m*(k-1);
        
        int lda = m*(k-1), nrhs = 1, ldb = m*(k-1), info;
        //Cholesky factorization of a symmetric (Hermitian) positive-definite matrix
        LAPACK_POTRF ( "U", &lda, XWX, &lda, &info );
        // Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite matrix
        LAPACK_POTRS ( "U", &lda, &nrhs, XWX, &lda, XWZ, &ldb, &info );  // Ab=solution vector
        if ( info != 0 )
            cout<<"error equation solver failed to converge: "<<info<<endl;
        for(int i=0;i<k-1;i++)
            for(int j=0;j<m;j++)
                x[i+j*(k-1)] = XWZ[j+i*m];
        
        //for(int i=0;i<m*(k-1);i++)
        //    cout<<x[i]<<" ";
        
        double sum0 = 0.0, sum1 = 0.0, normalizedDiff;
        for(int i=0;i<m*(k-1);i++)
        {
            sum0 += x[i] * x[i];
            sum1 += xOld[i] * xOld[i];
        }
        normalizedDiff = (sum0-sum1)/((double)m*(k-1));
        //cout<<" | "<<normalizedDiff<<endl;
        if(normalizedDiff < tolerance)  // break criteria
            return;
        
        for(int i=0;i<m*(k-1);i++)
            xOld[i] = x[i];
        
        for(int i=0;i<n*k;i++)
            etaOld[i] = eta[i];
        
        CBLAS_GEMM ( CblasRowMajor, CblasNoTrans, CblasNoTrans, n, k-1 , m, 1.0, A, m, x, k-1, 0.0, etaTmp, k-1 ); // etaTmp = A * x
        
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<k-1;j++)
                eta[j+i*k] = etaTmp[j+i*(k-1)];
            eta[k-1+i*k] = 0.0;
        }
        
        for(int backstep=0;backstep<11;backstep++)
        {
            for(int i=0;i<n*k;i++)
                pi[i] = exp(eta[i]);
            
            for(int i=0;i<n;i++)
            {
                REAL sum = 0.0;
                REAL* ptr = pi + i*k;
                for(int j=0;j<k;j++)
                    sum += ptr[j];
                for(int j=0;j<k;j++)
                    ptr[j] /= sum;
            }
            
            uint sum = 0;
            for(int i=0;i<n*k;i++)
                if(pi[i] > tolerance)
                    sum++;
            if(sum == n*k)
                break;
            else if(backstep < 10)
            {
                for(int i=0;i<n*k;i++)
                    eta[i] = etaOld[i] + (eta[i] - etaOld[i])*0.2;
            }
            else
            {
                for(int i=0;i<n*k;i++)
                    if(pi[i] <= tolerance)
                        pi[i] = tolerance;
                for(int i=0;i<n;i++)
                {
                    REAL sum = 0.0;
                    REAL* ptr = pi + i*k;
                    for(int j=0;j<k;j++)
                        sum += ptr[j];
                    for(int j=0;j<k;j++)
                        pi[j] /= sum;
                }
                for(int i=0;i<n*k;i++)
                    eta[i] = log(pi[i]);
                break;
            }
        }
    }
    
    delete[] xOld;
    delete[] pi;
    delete[] eta;
    delete[] etaTmp;
    delete[] etaOld;
    delete[] mu;
    delete[] XWX;
    delete[] XWZ;
    delete[] W;
    delete[] Z;
}

void NumericalTools::RidgeRegressionMultisolutionMulticall	(	REAL *	A,
		REAL *	b,
		REAL *	x,
		int	n,
		int	m,
		int	k,
		REAL	lambda,
		int	alloc = 0,
		bool	isRowMajor = false
	)			[static]

Ridge Regression with multiple target vectors (x precision)

Use this function, if you have to solve multiple linear equations.

A ... n x m matrix b ... n x k target vector x ... m x k solution vector

min ||A*x - b||

Definition at line 110 of file NumericalTools.cpp.

{
    static REAL *AA;
    static REAL *Ab;
    static REAL *AbTrans;
    static REAL *eigenvals;
    static REAL *work;
    static int *iwork;
    int oneInt=1, rank, info;
    REAL precision = 1e-15;
    static int lWork=-1;
    static REAL workSize[1];

    if ( alloc == ALLOC_MEM )
    {
        AA = new REAL[m*m];
        Ab = new REAL[m*k];
        AbTrans = new REAL[m*k];
        eigenvals = new REAL[m];
        lWork = -1;
        GELSS ( &m, &m, &k, AA, &m, Ab, &m, eigenvals, &precision, &rank, workSize, &lWork, &info );

        lWork = ( int ) ( workSize[0]+0.5 );
        work = new REAL[lWork];

        return;
    }

    if ( alloc == FREE_MEM )
    {
        if ( work )
            delete[] work;
        work = 0;
        if ( AA )
            delete[] AA;
        AA = 0;
        if ( Ab )
            delete[] Ab;
        Ab = 0;
        if ( AbTrans )
            delete[] AbTrans;
        AbTrans = 0;
        if ( eigenvals )
            delete[] eigenvals;
        eigenvals = 0;

        return;
    }

    if ( isRowMajor )
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, AA, m ); // AA = A'A
    else
        CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, AA, m ); // AA = A'A

    for ( int i=0; i < m; i++ )
        AA[i+i*m] += ( REAL ) n*lambda;  // A'A + lambda*I

    if ( isRowMajor )
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, k , n, 1.0, A, m, b, k, 0.0, Ab, k ); // Ab = A'b
    else
        CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b

    if ( isRowMajor )
    {
        for ( int i=0;i<k;i++ )
            for ( int j=0;j<m;j++ )
                AbTrans[j+i*m] = Ab[i+j*k];  // copy to colMajor
    }
    else
    {
        for ( int i=0;i<k;i++ )
            for ( int j=0;j<m;j++ )
                AbTrans[i+j*k] = Ab[i+j*k];  // copy to colMajor
    }

    GELSS ( &m, &m, &k, AA, &m, AbTrans, &m, eigenvals, &precision, &rank, work, &lWork, &info );

    if ( info != 0 )
        cout<<"error equation solver failed to converge: "<<info<<endl;

    if ( isRowMajor )
    {
        for ( int i=0;i<m;i++ ) // all features
            for ( int j=0;j<k;j++ ) // all classes (#outupts)
                x[j+i*k] = AbTrans[i+j*m];
    }
    else
    {
        for ( int i=0;i<m;i++ ) // all features
            for ( int j=0;j<k;j++ ) // all classes (#outupts)
                x[i+j*m] = AbTrans[i+j*m];
    }

}

void NumericalTools::RidgeRegressionMultisolutionSinglecall	(	REAL *	A,
		REAL *	b,
		REAL *	x,
		int	n,
		int	m,
		int	k,
		REAL	lambda,
		bool	isRowMajor,
		double *	rideModifier = 0
	)			[static]

Ridge Regression with multiple target vectors (x precision)

Use this function, if you have to solve multiple linear equations.

A ... n x m matrix b ... n x k target vector x ... m x k solution vector

min ||A*x - b||

Definition at line 218 of file NumericalTools.cpp.

{
    REAL *AA;
    REAL *Ab;
    REAL *AbTrans;
    REAL *eigenvals;
    //REAL *work;
    int *iwork;
    int oneInt=1, rank, info;
    REAL precision = 1e-15;
    int lWork=-1;
    REAL workSize[1];

    AA = new REAL[m*m];
    Ab = new REAL[m*k];
    AbTrans = new REAL[m*k];
    eigenvals = new REAL[m];
    lWork = -1;

    // get size of work array
    //GELSS(&m, &m, &k, AA, &m, Ab, &m, eigenvals, &precision, &rank, workSize, &lWork, &info);
    //lWork = (int)(workSize[0]+0.5);
    //work = new REAL[lWork];

    if ( isRowMajor )
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, AA, m ); // AA = A'A
    else
        CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, AA, m ); // AA = A'A

    if ( rideModifier )
        for ( int i=0; i < m; i++ )
            AA[i+i*m] += ( REAL ) n*lambda*rideModifier[i];  // A'A + lambda*I*modifier
    else
        for ( int i=0; i < m; i++ )
            AA[i+i*m] += ( REAL ) n*lambda;  // A'A + lambda*I

    if ( isRowMajor )
        CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, k , n, 1.0, A, m, b, k, 0.0, Ab, k ); // Ab = A'b
    else
        CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b

    if ( isRowMajor )
    {
        for ( int i=0;i<k;i++ )
            for ( int j=0;j<m;j++ )
                AbTrans[j+i*m] = Ab[i+j*k];  // copy to colMajor
    }
    else
    {
        for ( int i=0;i<k;i++ )
            for ( int j=0;j<m;j++ )
                AbTrans[i+j*k] = Ab[i+j*k];  // copy to colMajor
    }

    int lda = m, nrhs = k, ldb = m;

    //Cholesky factorization of a symmetric (Hermitian) positive-definite matrix
    LAPACK_POTRF ( "U", &m, AA, &lda, &info );

    // Solves a system of linear equations with a Cholesky-factored symmetric (Hermitian) positive-definite matrix
    LAPACK_POTRS ( "U", &m, &nrhs, AA, &lda, AbTrans, &ldb, &info );  // Y=solution vector


    // slower solver
    // Computes the minimum-norm solution to a linear least squares problem using the singular value decomposition of A
    //REAL *work;
    // get size of work array
    //GELSS(&m, &m, &k, AA, &m, Ab, &m, eigenvals, &precision, &rank, workSize, &lWork, &info);
    //lWork = (int)(workSize[0]+0.5);
    //work = new REAL[lWork];
    //GELSS(&m, &m, &k, AA, &m, AbTrans, &m, eigenvals, &precision, &rank, work, &lWork, &info);
    //delete[] work;

    if ( info != 0 )
        cout<<"error equation solver failed to converge: "<<info<<endl;

    if ( isRowMajor )
    {
        for ( int i=0;i<m;i++ ) // all features
            for ( int j=0;j<k;j++ ) // all classes (#outupts)
                x[j+i*k] = AbTrans[i+j*m];
    }
    else
    {
        for ( int i=0;i<m;i++ ) // all features
            for ( int j=0;j<k;j++ ) // all classes (#outupts)
                x[i+j*m] = AbTrans[i+j*m];
    }

    if ( AA )
        delete[] AA;
    AA = 0;
    if ( Ab )
        delete[] Ab;
    Ab = 0;
    if ( AbTrans )
        delete[] AbTrans;
    AbTrans = 0;
    if ( eigenvals )
        delete[] eigenvals;
    eigenvals = 0;
}

void NumericalTools::RidgeRegressionMultisolutionSinglecallGELSS	(	REAL *	A,
		REAL *	b,
		REAL *	x,
		int	n,
		int	m,
		int	k,
		REAL	lambda,
		bool	isRowMajor = false
	)			[static]

Ridge Regression with multiple target vectors (x precision) Solves the system with one GELSS call

Use this function, if you have to solve multiple linear equations.

A ... n x m matrix b ... n x k target vector x ... m x k solution vector

min ||A*x - b||

Definition at line 33 of file NumericalTools.cpp.

{
    int oneInt=1, rank, info;
    REAL precision = 1e-15;

    // the new one
    int nn = n + m;
    int lWork = -1;
    REAL workSize[1];
    REAL* AA = new REAL[nn*m];
    REAL* bb = new REAL[nn*k];
    REAL* eigenvals = new REAL[m];

    GELSS ( &nn, &m, &k, AA, &nn, bb, &nn, eigenvals, &precision, &rank, workSize, &lWork, &info );

    lWork = ( int ) ( workSize[0]+0.5 );
    REAL work[lWork];

    for ( int i=0;i<n;i++ )
    {
        for ( int j=0;j<m;j++ )
            AA[i + j*nn] = A[i + j*n];
        for ( int j=0;j<k;j++ )
            bb[i + j*nn] = b[i + j*n];
    }

    // add regularizer matrix to the bottom of A
    //
    // A = [1 2 3]  b = [1]
    //     [3 4 5]      [2]
    //     [.. ..]      [.]
    // lam*[1 0 0]      [0]  -> solving this extended equation system is
    // lam*[0 1 0]      [0]  -> equivalent to ridge regression !
    // lam*[0 0 1]      [0]  -> (Paper from Y.Koren 2007 used that solving schema)
    REAL lam = sqrt ( lambda* ( REAL ) n );
    for ( int i=0;i<m;i++ )
    {
        for ( int j=0;j<m;j++ )
        {
            if ( i == j )
                AA[j + i*nn + n] = lam;
            else
                AA[j + i*nn + n] = 0.0;
        }
        bb[i + n] = 0.0;
    }

    GELSS ( &nn, &m, &k, AA, &nn, bb, &nn, eigenvals, &precision, &rank, work, &lWork, &info );

    for ( int i=0;i<m;i++ )
        for ( int j=0;j<k;j++ )
            x[i+j*k] = bb[i+j*nn];

    if ( AA )
        delete[] AA;
    AA = 0;
    if ( bb )
        delete[] bb;
    bb = 0;
    if ( eigenvals )
        delete[] eigenvals;
    eigenvals = 0;
}

int NumericalTools::RidgeRegressionNonNegmulticall	(	REAL *	A,
		REAL *	b,
		REAL *	x,
		int	n,
		int	m,
		REAL	lamda,
		REAL	eps,
		int	maxIter,
		int	alloc = 0,
		bool	isRowMajor = true,
		int	debug = 0
	)			[static]

Solves the system Ax=b with x>0 The n x m matrix A is stored row-wise. So A[1] is A_01 Call the first time with alloc=ALLOC_MEM, then the internal static fields are allocated All arrays REAL single precision

Parameters:

	A	The n x m matrix (n rows, m columns)
	b	The target vector of length n
	x	The solution vector of length m
	n	The dimension n (rows)
	m	The dimension m (cols)
	eps	Precision of the norm from r (e.g. 1e-6 .. 1e-10)
	maxIter	Max. number of solver iteration
	alloc	Two cases, 1) ALLOC_MEM: allocates internal field and returns 2) FREE_MEM frees the memory and returns
	debug	If set to 1 then the fields are printed to cout

Returns:

The number of iterations for convergence

Definition at line 700 of file NumericalTools.cpp.

{
    static REAL *r, *tmp, *ATA, *Ab;

    if ( alloc == ALLOC_MEM )
    {
        cout<<"Allocate memory in nonNeg solver"<<endl;
        r = new REAL[m];
        tmp = new REAL[m];
        ATA = new REAL[m*m];
        Ab = new REAL[m];
        return -1;
    }
    if ( alloc == FREE_MEM )
    {
        cout<<"Free memory in nonNeg solver"<<endl;
        if ( r )
            delete[] r;
        r = 0;
        if ( tmp )
            delete[] tmp;
        tmp = 0;
        if ( ATA )
            delete[] ATA;
        ATA = 0;
        if ( Ab )
            delete[] Ab;
        Ab = 0;
        return -1;
    }

    int nIter = 0;
    unsigned int i;
    REAL alpha = 0.0, tmp0, tmp1;
    REAL norm = 1.0;

    if ( lamda > 0.0 )
    {
        if ( isRowMajor )
            CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, ATA, m ); // AA = A'A
        else
            CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, ATA, m ); // AA = A'A

        if ( isRowMajor )
            CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, m, b, 1, 0.0, Ab, 1 ); // calc Ab = A'*B
        else
            CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b

        // add constant diagonal elements for regularization
        for ( i=0;i<m;i++ )
            ATA[i+i*m] += lamda * ( REAL ) n;

    }
    else
    {
        memcpy ( ( char* ) ATA, ( char* ) A, sizeof ( REAL ) *m*m );
        memcpy ( ( char* ) Ab, ( char* ) b, sizeof ( REAL ) *m );
    }

    // init solution
    for ( i=0;i<m;i++ )
        x[i] = r[i] = 1.0;

    while ( norm>eps )
    {
        CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, x, 1, 0.0, tmp, 1 );  // tmp = ATA*x
        for ( i=0;i<m;i++ )
            r[i] = Ab[i] - tmp[i];
        for ( i=0;i<m;i++ )
            if ( x[i]==0.0 &&  r[i]<0.0 )
                r[i] = 0.0;
        CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, r, 1, 0.0, tmp, 1 );  // tmp = A*r
        tmp0 = tmp1 = 0.0;
        for ( i=0;i<m;i++ )
        {
            tmp0 += r[i]*r[i];   // nominator
            tmp1 += r[i]*tmp[i]; // denominator
        }
        if ( tmp0!=0.0 && tmp1!=0.0 )
            alpha = tmp0 / tmp1;  // alpha = r'r / r'Ar   (tmp=Ar)
        for ( i=0;i<m;i++ )
        {
            if ( r[i]<0.0 )
            {
                tmp0 = -x[i]/r[i];
                alpha = alpha<tmp0? alpha : tmp0;
            }
        }
        for ( i=0;i<m;i++ )
            x[i] += alpha * r[i];
        norm = 0.0;
        for ( i=0;i<m;i++ )
        {
            if ( x[i]<1e-9 )
                x[i] = 0.0;  // x(x<1e-6) = 0
            norm += r[i]*r[i];
        }
        norm = sqrt ( norm );
        nIter++;

        // break criteria, if no convergence
        if ( nIter>=maxIter )
            break;

        if ( debug )
        {
            cout<<"norm="<<norm<<"  alpha="<<alpha<<endl;
            cout<<"r=";
            for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
            cout<<endl;
            cout<<"x=";
            for ( i=0;i<m;i++ ) cout<<x[i]<<" ";
            cout<<endl;
            cout<<"tmp=";
            for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
            cout<<endl;
            cout<<endl;
        }
    }

    // check the solution for nan of inf values
    for ( i=0;i<m;i++ )
        if ( isnan ( x[i] ) || isinf ( x[i] ) )
            cout<<"DIVERGENCE of nonNegSolver (nIter="<<nIter<<")"<<endl;

    return nIter;
}

int NumericalTools::RidgeRegressionNonNegSinglecall	(	REAL *	A,
		REAL *	b,
		REAL *	x,
		int	n,
		int	m,
		REAL	lamda,
		REAL	eps,
		int	maxIter,
		bool	isRowMajor = true,
		int	debug = 0
	)

Solves the system Ax=b with x>0 The n x m matrix A is stored row-wise. So A[1] is A_01 Call the first time with alloc=ALLOC_MEM, then the internal static fields are allocated All arrays REAL single precision

Parameters:

	A	The n x m matrix (n rows, m columns)
	b	The target vector of length n
	x	The solution vector of length m
	n	The dimension n (rows)
	m	The dimension m (cols)
	eps	Precision of the norm from r (e.g. 1e-6 .. 1e-10)
	maxIter	Max. number of solver iteration
	debug	If set to 1 then the fields are printed to cout

Returns:

The number of iterations for convergence

Definition at line 562 of file NumericalTools.cpp.

{
    static REAL *r, *tmp, *ATA, *Ab;

    r = new REAL[m];
    tmp = new REAL[m];
    ATA = new REAL[m*m];
    Ab = new REAL[m];

    int nIter = 0;
    unsigned int i;
    REAL alpha = 0.0, tmp0, tmp1;
    REAL norm = 1.0;

    if ( lamda > 0.0 )
    {
        if ( isRowMajor )
            CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, m, A, m, 0.0, ATA, m ); // AA = A'A
        else
            CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, m, n, 1.0, A, n, A, n, 0.0, ATA, m ); // AA = A'A

        if ( isRowMajor )
            CBLAS_GEMM ( CblasRowMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, m, b, 1, 0.0, Ab, 1 ); // calc Ab = A'*B
        else
            CBLAS_GEMM ( CblasColMajor, CblasTrans, CblasNoTrans, m, 1 , n, 1.0, A, n, b, n, 0.0, Ab, m ); // Ab = A'b

        // add constant diagonal elements for regularization
        for ( i=0;i<m;i++ )
            ATA[i+i*m] += lamda * ( REAL ) n;

    }
    else
    {
        memcpy ( ( char* ) ATA, ( char* ) A, sizeof ( REAL ) *m*m );
        memcpy ( ( char* ) Ab, ( char* ) b, sizeof ( REAL ) *m );
    }

    // init solution
    for ( i=0;i<m;i++ )
        x[i] = r[i] = 1.0;

    while ( norm>eps )
    {
        CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, x, 1, 0.0, tmp, 1 );  // tmp = ATA*x
        for ( i=0;i<m;i++ )
            r[i] = Ab[i] - tmp[i];
        for ( i=0;i<m;i++ )
            if ( x[i]==0.0 &&  r[i]<0.0 )
                r[i] = 0.0;
        CBLAS_GEMV ( CblasRowMajor, CblasNoTrans, m, m, 1.0, ATA, m, r, 1, 0.0, tmp, 1 );  // tmp = A*r
        tmp0 = tmp1 = 0.0;
        for ( i=0;i<m;i++ )
        {
            tmp0 += r[i]*r[i];   // nominator
            tmp1 += r[i]*tmp[i]; // denominator
        }
        if ( tmp0!=0.0 && tmp1!=0.0 )
            alpha = tmp0 / tmp1;  // alpha = r'r / r'Ar   (tmp=Ar)
        for ( i=0;i<m;i++ )
        {
            if ( r[i]<0.0 )
            {
                tmp0 = -x[i]/r[i];
                alpha = alpha<tmp0? alpha : tmp0;
            }
        }
        for ( i=0;i<m;i++ )
            x[i] += alpha * r[i];
        norm = 0.0;
        for ( i=0;i<m;i++ )
        {
            if ( x[i]<1e-9 )
                x[i] = 0.0;  // x(x<1e-6) = 0
            norm += r[i]*r[i];
        }
        norm = sqrt ( norm );
        nIter++;

        // break criteria, if no convergence
        if ( nIter>=maxIter )
            break;

        if ( debug )
        {
            cout<<"norm="<<norm<<"  alpha="<<alpha<<endl;
            cout<<"r=";
            for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
            cout<<endl;
            cout<<"x=";
            for ( i=0;i<m;i++ ) cout<<x[i]<<" ";
            cout<<endl;
            cout<<"tmp=";
            for ( i=0;i<m;i++ ) cout<<r[i]<<" ";
            cout<<endl;
            cout<<endl;
        }
    }

    // check the solution for nan of inf values
    for ( i=0;i<m;i++ )
        if ( isnan ( x[i] ) || isinf ( x[i] ) )
            cout<<"DIVERGENCE of nonNegSolver (nIter="<<nIter<<")"<<endl;

    if ( r )
        delete[] r;
    r = 0;
    if ( tmp )
        delete[] tmp;
    tmp = 0;
    if ( ATA )
        delete[] ATA;
    ATA = 0;
    if ( Ab )
        delete[] Ab;
    Ab = 0;

    return nIter;
}

---

The documentation for this class was generated from the following files:

• ELF/NumericalTools.h
• ELF/NumericalTools.cpp