Numerical Tools

Collaboration diagram for Numerical Tools:

Classes
class	GaussianInterpolator

Functions
void	randomPermutation (int N, MultidimArray< int > &result)
void	powellOptimizer (Matrix1D< double > &p, int i0, int n, double(*f)(double *, void *), void *prm, double ftol, double &fret, int &iter, const Matrix1D< double > &steps, bool show=false)
double	solveNonNegative (const Matrix2D< double > &A, const Matrix1D< double > &b, Matrix1D< double > &result)
void	solveViaCholesky (const Matrix2D< double > &A, const Matrix1D< double > &b, Matrix1D< double > &result)
void	evaluateQuadratic (const Matrix1D< double > &x, const Matrix1D< double > &c, const Matrix2D< double > &H, double &val, Matrix1D< double > &grad)
void	quadraticProgramming (const Matrix2D< double > &C, const Matrix1D< double > &d, const Matrix2D< double > &A, const Matrix1D< double > &b, const Matrix2D< double > &Aeq, const Matrix1D< double > &beq, Matrix1D< double > &bl, Matrix1D< double > &bu, Matrix1D< double > &x)
void	leastSquare (const Matrix2D< double > &C, const Matrix1D< double > &d, const Matrix2D< double > &A, const Matrix1D< double > &b, const Matrix2D< double > &Aeq, const Matrix1D< double > &beq, Matrix1D< double > &bl, Matrix1D< double > &bu, Matrix1D< double > &x)
void	regularizedLeastSquare (const Matrix2D< double > &A, const Matrix1D< double > &d, double lambda, const Matrix2D< double > &G, Matrix1D< double > &x)
template<typename T >
void	solve (const Matrix2D< T > &A, const Matrix1D< T > &b, Matrix1D< T > &result)
template<typename T >
void	solveBySVD (const Matrix2D< T > &A, const Matrix1D< T > &b, Matrix1D< double > &result, double tolerance)
template<typename T >
void	solve (const Matrix2D< T > &A, const Matrix2D< T > &b, Matrix2D< T > &result)

Detailed Description

Function Documentation

evaluateQuadratic()

void evaluateQuadratic	(	const Matrix1D< double > &	x,
		const Matrix1D< double > &	c,
		const Matrix2D< double > &	H,
		double &	val,
		Matrix1D< double > &	grad
	)

Evaluate quadratic form

Given x, c and H this function returns the value of the quadratic form val=c^t*x+0.5*x^t*H^t*H*x and the gradient of the quadratic form at x grad=c+H*x.

Exceptions are thrown if the vectors and matrices do not have consistent dimensions.

Definition at line 126 of file numerical_tools.cpp.

{
    if (x.size() != c.size())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Eval_quadratic: Not compatible sizes in x and c");
    if (H.Xdim() != x.size())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Eval_quadratic: Not compatible sizes in x and H");

    // H*x, store in grad
    grad.initZeros(x.size());
    for (size_t i = 0; i < H.Ydim(); i++)
        for (size_t j = 0; j < x.size(); j++)
            grad(i) += H(i, j) * x(j);

    // Now, compute c^t*x+1/2*x^t*H*x
    // Add c to the gradient
    double quad = 0;
    val = 0;
    for (size_t j = 0; j < x.size(); j++)
    {
        quad += grad(j) * grad(j); // quad=x^t*H^t*H*x
        val += c(j) * x(j);  // val=c^t*x

        grad(j) += c(j);     // grad+=c
    }
    val += 0.5 * quad;
}

leastSquare()

void leastSquare	(	const Matrix2D< double > &	C,
		const Matrix1D< double > &	d,
		const Matrix2D< double > &	A,
		const Matrix1D< double > &	b,
		const Matrix2D< double > &	Aeq,
		const Matrix1D< double > &	beq,
		Matrix1D< double > &	bl,
		Matrix1D< double > &	bu,
		Matrix1D< double > &	x
	)

Solves the least square problem

min 0.5*(Norm(C*x-d))   subject to:  A*x <= b
x                                    Aeq*x=beq
                                     bl<=x<=bu

Definition at line 322 of file numerical_tools.cpp.

{
    // Convert d to Matrix2D for multiplication
    Matrix2D<double> P;
    P.fromVector(d);
    P = -2.0 * (P.transpose() * C);
    P = P.transpose();

    //Convert back to vector for passing it to quadraticProgramming
    Matrix1D<double> newd;
    P.toVector(newd);

    quadraticProgramming(C.transpose()*C, newd, A, b, Aeq, beq, bl, bu, x);
}

powellOptimizer()

void powellOptimizer	(	Matrix1D< double > &	p,
		int	i0,
		int	n,
		double(*)(double *, void *)	f,
		void *	prm,
		double	ftol,
		double &	fret,
		int &	iter,
		const Matrix1D< double > &	steps,
		bool	show = false
	)

Optimize using Powell's method.

See Numerical Recipes Chapter 10.

Problem: Minimize f(x) starting at point p. n is the dimension of x. If changes are smaller than ftol then stop. The number of iterations is returned in iter, fret contains the function value at the minimum and p contains the minimum.

Watch out that the evaluating function f must consider that x starts at index 1, at least, and goes until n. i0 is used to allow optimization in large spaces where only one part is going to be optimized. Thus, if in a vector of dimension 20 you want to optimize the first 3 components then i0=1, n=3; if you want to optimize it all, i0=1, n=20; and if you want to optimize the last five components i0=15, n=5.

The steps define the allowed steps on each variable. When a variable has large ranges this step should be larger to allow bigger movements. The steps should have only the dimension of the optimized part (3,20 and 5 in the examples above).

The option show forces the routine to show the convergence path

If your function needs extra parameters you can pass them through the void pointer prm. If you don't need this feature set it to NULL.

Example of use:

MultidimArray<double> x(8), steps(8);
double fitness;
int iter;
steps.initConstant(1);
x.initZeros();
powellOptimizer(x,1,8,&wrapperFitness,NULL,0.01,fitness,iter,steps,true);

quadraticProgramming()

void quadraticProgramming	(	const Matrix2D< double > &	C,
		const Matrix1D< double > &	d,
		const Matrix2D< double > &	A,
		const Matrix1D< double > &	b,
		const Matrix2D< double > &	Aeq,
		const Matrix1D< double > &	beq,
		Matrix1D< double > &	bl,
		Matrix1D< double > &	bu,
		Matrix1D< double > &	x
	)

Solves Quadratic programming subproblem

min 0.5*x'Cx + d'x   subject to:  A*x <= b
 x                                Aeq*x=beq
                                  bl<=x<=bu

Definition at line 222 of file numerical_tools.cpp.

{
    CDAB prm;
    prm.C = C;
    prm.D.fromVector(d);
    prm.A.initZeros(A.Ydim() + Aeq.Ydim(), A.Xdim());
    prm.B.initZeros(prm.A.Ydim(), 1);


    // Copy Inequalities
    for (size_t i = 0; i < A.Ydim(); i++)
    {
        for (size_t j = 0; j < A.Xdim(); j++)
            prm.A(i, j) = A(i, j);
        prm.B(i, 0) = b(i);
    }

    // Copy Equalities
    for (size_t i = 0; i < Aeq.Ydim(); i++)
    {
        for (size_t j = 0; j < Aeq.Xdim(); j++)
            prm.A(i + A.Ydim(), j) = Aeq(i, j);
        prm.B(i + A.Ydim(), 0) = beq(i);
    }

    double bigbnd = 1e30;
    // Bounds
    if (bl.size() == 0)
    {
        bl.resize(C.Xdim());
        bl.initConstant(-bigbnd);
    }
    if (bu.size() == 0)
    {
        bu.resize(C.Xdim());
        bu.initConstant(bigbnd);
    }

    // Define intermediate variables
    int    mode = 100;  // CFSQP mode
    int    iprint = 0;  // Debugging
    int    miter = 1000;  // Maximum number of iterations
    double eps = 1e-4; // Epsilon
    double epsneq = 1e-4; // Epsilon for equalities
    double udelta = 0.e0; // Finite difference approximation
    // of the gradients. Not used in this function
    int    nparam = C.Xdim(); // Number of variables
    int    nf = 1;          // Number of objective functions
    int    neqn = Aeq.Ydim();        // Number of nonlinear equations
    int    nineqn = A.Ydim();      // Number of nonlinear inequations
    int    nineq = A.Ydim();  // Number of linear inequations
    int    neq = Aeq.Ydim();  // Number of linear equations
    int    inform;
    int    ncsrl = 0;
    int    ncsrn = 0;
    int    nfsr = 0;
    int    mesh_pts[] = {0};

    if (x.size() == 0)
        x.initZeros(nparam);
    Matrix1D<double> f(nf);
    Matrix1D<double> g(nineq + neq);
    Matrix1D<double> lambda(nineq + neq + nf + nparam);

    // Call the minimization routine
    cfsqp(nparam, nf, nfsr, nineqn, nineq, neqn, neq, ncsrl, ncsrn, mesh_pts,
          mode, iprint, miter, &inform, bigbnd, eps, epsneq, udelta,
          MATRIX1D_ARRAY(bl), MATRIX1D_ARRAY(bu),
          MATRIX1D_ARRAY(x),
          MATRIX1D_ARRAY(f), MATRIX1D_ARRAY(g),
          MATRIX1D_ARRAY(lambda),
          //  quadprg_obj32,quadprog_cntr32,quadprog_grob32,quadprog_grcn32,
          quadraticProgramming_obj32, quadraticProgramming_cntr32, grobfd, grcnfd,
          (void*)&prm);

#ifdef DEBUG

    if (inform == 0)
        std::cout << "SUCCESSFUL RETURN. \n";
    if (inform == 1 || inform == 2)
        std::cout << "\nINITIAL GUESS INFEASIBLE.\n";
    if (inform == 3)
        printf("\n MAXIMUM NUMBER OF ITERATIONS REACHED.\n");
    if (inform > 3)
        printf("\ninform=%d\n", inform);
#endif
}

randomPermutation()

void randomPermutation	(	int	N,
		MultidimArray< int > &	result
	)

Generate random permutation

Generate a random permutation of the numbers between 0 and N-1

Definition at line 33 of file numerical_tools.cpp.

{
    MultidimArray<double> aux;
    aux.resize(N);
    aux.initRandom(0,1);

    aux.indexSort(result);
    result-=1;
}

regularizedLeastSquare()

void regularizedLeastSquare	(	const Matrix2D< double > &	A,
		const Matrix1D< double > &	d,
		double	lambda,
		const Matrix2D< double > &	G,
		Matrix1D< double > &	x
	)

Solves the regularized least squares problem

min Norm(A*x-d) + lambda * norm (G*x)

Give an empty G matrix (NULL matrix) if G is the identity matrix If AtA is not a NULL matrix, then the product AtA is not computed.

Definition at line 342 of file numerical_tools.cpp.

{
    int Nx=A.Xdim(); // Number of variables

    Matrix2D<double> X(Nx,Nx); // X=(A^t * A +lambda *G^t G)
    // Compute A^t*A
    FOR_ALL_ELEMENTS_IN_MATRIX2D(X)
    // Compute the dot product of the i-th and j-th columns of A
    for (size_t k=0; k<A.Ydim(); k++)
        X(i,j) += A(k,i) * A(k,j);

    // Compute lambda*G^t*G
    if (G.Xdim()==0)
        for (int i=0; i<Nx; i++)
            X(i,i)+=lambda;
    else
        FOR_ALL_ELEMENTS_IN_MATRIX2D(X)
        // Compute the dot product of the i-th and j-th columns of G
        for (size_t k=0; k<G.Ydim(); k++)
            X(i,j) += G(k,i) * G(k,j);

    // Compute A^t*d
    Matrix1D<double> Atd(Nx);
    FOR_ALL_ELEMENTS_IN_MATRIX1D(Atd)
    // Compute the dot product of the i-th column of A and d
    for (size_t k=0; k<A.Ydim(); k++)
        Atd(i) += A(k,i) * d(k);

    // Compute the inverse of X
    Matrix2D<double> Xinv;
    X.inv(Xinv);

    // Now multiply Xinv * A^t * d
    matrixOperation_Ax(Xinv, Atd, x);
}

solve() [1/2]

template<typename T >

void solve	(	const Matrix2D< T > &	A,
		const Matrix1D< T > &	b,
		Matrix1D< T > &	result
	)

Solves a linear equation system by LU decomposition.

Definition at line 215 of file numerical_tools.h.

{
    if (MAT_XSIZE(A) == 0)
        REPORT_ERROR(ERR_MATRIX_EMPTY, "Solve: Matrix is empty");

    if (MAT_XSIZE(A) != MAT_YSIZE(A))
        REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Matrix is not squared");

    if (MAT_XSIZE(A) != b.size())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Different sizes of Matrix and Vector");

    if (b.isRow())
        REPORT_ERROR(ERR_MATRIX_DIM, "Solve: Not correct vector shape");

    // Perform LU decomposition and then solve
    Matrix1D< int > indx;
    T d;
    Matrix2D<T> LU;
    ludcmp(A, LU, indx, d);
    result = b;
    lubksb(LU, indx, result);
}

solve() [2/2]

template<typename T >

void solve	(	const Matrix2D< T > &	A,
		const Matrix2D< T > &	b,
		Matrix2D< T > &	result
	)

Solves a linear equation system by Gaussian elimination.

Definition at line 2218 of file numerical_tools.cpp.

{
    if (A.Xdim() == 0)
        REPORT_ERROR(ERR_MATRIX_EMPTY, "Solve: Matrix is empty");

    if (A.Xdim() != A.Ydim())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Matrix is not squared");

    if (A.Ydim() != b.Ydim())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Different sizes of A and b");

    // Solve
    result = b;
    Matrix2D<T> Aux = A;
    gaussj(Aux.adaptForNumericalRecipes2(), Aux.Ydim(),
           result.adaptForNumericalRecipes2(), b.Xdim());
}

solveBySVD()

template<typename T >

void solveBySVD	(	const Matrix2D< T > &	A,
		const Matrix1D< T > &	b,
		Matrix1D< double > &	result,
		double	tolerance
	)

Solves a linear equation system by SVD decomposition.

Definition at line 2179 of file numerical_tools.cpp.

{
    if (A.Xdim() == 0)
        REPORT_ERROR(ERR_MATRIX_EMPTY, "Solve: Matrix is empty");

    if (A.Xdim() != A.Ydim())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Matrix is not squared");

    if (A.Xdim() != b.size())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Solve: Different sizes of Matrix and Vector");

    if (b.isRow())
        REPORT_ERROR(ERR_MATRIX_DIM, "Solve: Not correct vector shape");

    // First perform de single value decomposition
    // Xmipp interface that calls to svdcmp of numerical recipes
    Matrix2D< double > u;
    Matrix2D< double > v;
    Matrix1D< double > w;
    svdcmp(A, u, w, v);

    // Here is checked if eigenvalues of the svd decomposition are acceptable
    // If a value is lower than tolerance, the it's zeroed, as this increases
    // the precision of the routine.
    for (int i = 0; i < w.size(); i++)
        if (w(i) < tolerance)
            w(i) = 0;

    // Set size of matrices
    result.resize(b.size());

    // Xmipp interface that calls to svdksb of numerical recipes
    Matrix1D< double > bd;
    typeCast(b, bd);
    svbksb(u, w, v, bd, result);
}

solveNonNegative()

double solveNonNegative	(	const Matrix2D< double > &	A,
		const Matrix1D< double > &	b,
		Matrix1D< double > &	result
	)

Solve equation system, nonnegative solution

The equation system is defined by Ax=b, it is solved for x. x is forced to be nonnegative. It is designed to cope with large equation systems. This function is borrowed from LAPACK nnls.

The norm of the vector Ax-b is returned.

Definition at line 82 of file numerical_tools.cpp.

{
    if (C.Xdim() == 0)
        REPORT_ERROR(ERR_MATRIX_EMPTY, "Solve_nonneg: Matrix is empty");
    if (C.Ydim() != d.size())
        REPORT_ERROR(ERR_MATRIX_SIZE, "Solve_nonneg: Different sizes of Matrix and Vector");
    if (d.isRow())
        REPORT_ERROR(ERR_MATRIX_DIM, "Solve_nonneg: Not correct vector shape");

    Matrix2D<double> Ct;
    Ct=C.transpose();

    result.initZeros(Ct.Ydim());
    double rnorm;

    // Watch out that matrix Ct is transformed.
    int success = nnls(MATRIX2D_ARRAY(Ct), Ct.Xdim(), Ct.Ydim(),
                       MATRIX1D_ARRAY(d),
                       MATRIX1D_ARRAY(result),
                       &rnorm, nullptr, nullptr, nullptr);
    if (success == 1)
        std::cerr << "Warning, too many iterations in nnls\n";
    else if (success == 2)
        REPORT_ERROR(ERR_MEM_NOTENOUGH, "Solve_nonneg: Not enough memory");
    return rnorm;
}

solveViaCholesky()

void solveViaCholesky	(	const Matrix2D< double > &	A,
		const Matrix1D< double > &	b,
		Matrix1D< double > &	result
	)

Solve equation system, symmetric positive-definite matrix

The equation system is defined by Ax=b, it is solved for x. This method can only be applied if A is positive-definite matrix and symmetric. It applies a Cholesky factorization and backsubstitution (see Numerical Recipes).

Definition at line 111 of file numerical_tools.cpp.

{
    Matrix2D<double> Ap = A;
    Matrix1D<double> p(A.Xdim());
    result.resize(A.Xdim());
    choldc(Ap.adaptForNumericalRecipes2(), A.Xdim(),
           p.adaptForNumericalRecipes());
    cholsl(Ap.adaptForNumericalRecipes2(), A.Xdim(),
           p.adaptForNumericalRecipes(), b.adaptForNumericalRecipes(),
           result.adaptForNumericalRecipes());
}