kSVD

Collaboration diagram for kSVD:

Functions
double	orthogonalMatchingPursuit (const Matrix1D< double > &x, const Matrix2D< double > &D, int S, Matrix1D< double > &alpha)
double	lasso (const Matrix1D< double > &x, const Matrix2D< double > &D, const Matrix2D< double > &DtD, const Matrix2D< double > &DtDlambdaInv, double lambda, Matrix1D< double > &alpha, const int maxIter=20, const double tol=0.005)

Variables
constexpr int	OMP_PROJECTION = 1
constexpr int	LASSO_PROJECTION = 2

Detailed Description

Function Documentation

lasso()

double lasso	(	const Matrix1D< double > &	x,
		const Matrix2D< double > &	D,
		const Matrix2D< double > &	DtD,
		const Matrix2D< double > &	DtDlambdaInv,
		double	lambda,
		Matrix1D< double > &	alpha,
		const int	maxIter = 20,
		const double	tol = 0.005
	)

Lasso projection

This function looks for the best approximation with as few components as possible (regularized through a L1 norm) of a given dictionary D. The representation is alpha, i.e., x=D*alpha.

DtD is D^t*D and must be precomputed before calling the dictionary

The approximation error is returned.

This implementation is based on Mark Schmidt http://www.cs.ubc.ca/~schmidtm/Software/lasso.html

Definition at line 149 of file kSVD.cpp.

{
    // Compute the ridge least squares solution
    // Compute D^t*x
    Matrix1D<double> Dtx;
    matrixOperation_Atx(D,x,Dtx); // Dtx=D^t*x
    matrixOperation_Ax(DtDlambdaInv,Dtx,alpha); // alpha=DtDlambdaInv*Dtx

    // Iterate over alpha
    Matrix1D<double> alphaOld;
    int iter=0;
    bool finished=false;
    while (!finished)
    {
        // Keep the current alpha
        alphaOld=alpha;
        
        // Update alpha
        FOR_ALL_ELEMENTS_IN_MATRIX1D(alpha)
        {
            double S=-Dtx(i);
            for (size_t j=0; j<alpha.size(); j++)
                if (i!=j) S+=DtD(i,j)*alpha(j);
            if (fabs(S)<lambda)
                alpha(i)=0;
            else if (S>lambda)
                alpha(i)=(lambda-S)/DtD(i,i);
            else
                alpha(i)=(-lambda-S)/DtD(i,i);
        }
        
        // Prepare for next iteration
        iter++;
        double normDiff=0;
        double normAlpha=0;
        FOR_ALL_ELEMENTS_IN_MATRIX1D(alpha)
        {
            double diff=alpha(i)-alphaOld(i);
            normDiff+=diff*diff;
            normAlpha+=alpha(i)*alpha(i);
        }
        finished=(iter>=maxIter) || normDiff/normAlpha<tol;
    }
    
    // Compute the approximation error
    Matrix1D<double> xp(x.size());
    for (size_t j=0; j<D.Xdim(); j++)
        if (alpha(j)!=0)
            for (size_t i=0; i<D.Ydim(); i++)
                xp(i)+=D(i,j)*alpha(j);
    double approximationError=0;
    FOR_ALL_ELEMENTS_IN_MATRIX1D(xp)
    {
        double diff=xp(i)-x(i);
        approximationError+=diff*diff;
    }
    
    return sqrt(ABS(approximationError));
}

orthogonalMatchingPursuit()

double orthogonalMatchingPursuit	(	const Matrix1D< double > &	x,
		const Matrix2D< double > &	D,
		int	S,
		Matrix1D< double > &	alpha
	)

Orthogonal matching pursuit (OMP)

This function looks for the best approximation with only S components of a given dictionary D. The representation is alpha, i.e., x=D*alpha.

The approximation error is returned.

This implementation is based on Karl Skretting http://www.ux.uis.no/~karlsk/proj02/index.html

Definition at line 31 of file kSVD.cpp.

{
    // Compute approximation error
    double approximationError=0;
    FOR_ALL_ELEMENTS_IN_MATRIX1D(x)
        approximationError+=x(i)*x(i);
    double errorLimit=approximationError*1e-10;

    // Initialize the number of available atoms
    int N=D.Ydim(); // Dimension of each atom
    if (0 == N) REPORT_ERROR(ERR_PARAM_INCORRECT, "Ydim of each atom must not be zero (0)");
    int K=D.Xdim(); // Number of atoms in the dictionary
    Matrix1D<int> availableAtoms;
    availableAtoms.initZeros(K);
    availableAtoms.initConstant(1);
    
    // Compute the projection of x onto each of the atoms and look for the
    // maximum
    Matrix1D<double> c;
    Matrix1D<double> w;
    Matrix1D<double> e;
    Matrix1D<double> u;
    c.initZeros(K);
    e.initZeros(K); e.initConstant(1);
    u.initZeros(K); u.initConstant(1);

    // Compute dot product between x and di
    int km=0;
    double ckm=-1;
    for (int k=0; k<K; k++)
    {
        for (int j=0; j<N; j++)
            c(k)+=D(j,k)*x(j);
        
        // Check if this is the maximum
        if (ABS(c(k))>ckm)
        {
            km=k;
            ckm=ABS(c(k));
        }
    }
    
    // Delete Js from the available indexes and update the approximation error
    availableAtoms(km)=0;
    approximationError-=ckm*ckm;
    
    // Annotate Jk
    std::vector<int> J;
    int s=0;
    J.push_back(km);

    // Build the R matrix
    Matrix2D<double> R(S,K);
    R(0,km)=u(km);

    while (s<S-1 && approximationError>errorLimit)
    {
        // Update all the rest of R except Js
        ckm=0;
        int nextKm=-1;
        for (int k=0; k<K; k++)
        {
            // If already used, skip it
            if (!availableAtoms(k)) continue;
        
            // Compute the product between the atoms km and k
            double aux=0;
            for (int j=0; j<N; j++)
                aux+=D(j,km)*D(j,k);
            R(s,k)=aux;
            
            // Readjust R
            for (int n=0; n<s; n++)
                R(s,k)-=R(n,km)*R(n,k);
            if (u(km)>XMIPP_EQUAL_ACCURACY)
                R(s,k)/=u(km);
            
            // Update the rest of variables
            c(k)=c(k)*u(k)-c(km)*R(s,k);
            e(k)-=R(s,k)*R(s,k);
            u(k)=sqrt(ABS(e(k)));
            if (u(k)!=0) c(k)/=u(k);

            // Check if this is the best c(k)
            if (ABS(c(k))>ckm)
            {
                nextKm=k;
                ckm=ABS(c(k));
            }
        }
        
        // Delete km from the available indexes and update the approximation error
        km=nextKm;
        J.push_back(km);
        availableAtoms(km)=0;
        s++;
        R(s,km)=u(km);
        approximationError-=ckm*ckm;
    }

    // Perform backsubstitution
    alpha.initZeros(K);
    for (int k=s; k>=0; k--)
    {
        int Jk=J[k];
        for (int n=s; n>=k+1; n--)
            c(Jk)-=R(k,J[n])*c(J[n]);
        if (R(k,Jk)!=0) c(Jk)/=R(k,Jk);
        alpha(Jk)=c(Jk);
    }
    
    return sqrt(ABS(approximationError));
}

Variable Documentation

LASSO_PROJECTION

constexpr int LASSO_PROJECTION = 2

Definition at line 72 of file kSVD.h.

OMP_PROJECTION

constexpr int OMP_PROJECTION = 1

Definition at line 71 of file kSVD.h.