pls.cpp

Go to the documentation of this file.

/* *********************
 * @file pls.cpp
 *
 * @author Stefan Ulbrich
 * @date 2013-2014
 *
 * @brief This file contains implementation of the PLS algorithm.
 * *********************/
#include "kbm.h"
#include <vector>
 
 
namespace memoryx::KBM
{
    /*
     * In case that regression and reconstruction should be separated later,
     * store the projection parameter in the data structure defined if
     * the makro is activated. Will come in handy for an incremental version.
     */
    //#define SEPARATE_REGRESSION_RECONSTRUCTION
 
#ifdef SEPARATE_REGRESSION_RECONSTRUCTION
    struct Projection
    {
        Vector p;
        Vector u;
        Real q;
        Projection(Vector p_, Vector u_, Real q_) : p(p_), u(u_), q(q_) {};
    };
#endif
 
    Matrix PLS::solve(const Matrix& input, const Matrix output, Real threshold)
    {
 
#ifdef SEPARATE_REGRESSION_RECONSTRUCTION
        std::vector<std::vector<Projection> > projections;
#endif
 
        // Results stored in the following matrix. It has to be initialized with zeros.
        Matrix parameters = Matrix::Zero(output.rows(), input.rows());
 
        // PLS 1 optimizes all output dimensions separately
        for (int d = 0; d < output.rows(); d++)
        {
 
#ifdef SEPARATE_REGRESSION_RECONSTRUCTION
            projections.push_back(std::vector<Projection>());
#endif
            // Make working copies
            Matrix input_d = input;
            Vector output_d = output.row(d);
 
            // In order to reconstruct the parameters, the identity matrix is transformed according to
            // projections defined during the regression.
            Matrix identity = Matrix::Identity(input.rows(), input.rows());
 
            for (int k = 0; k < input.rows(); k++) // maximal 3^{d_i} projections
            {
                // Step 1: Regression
                Vector u = input_d * output_d;
                Vector s = input_d.transpose() * u;
                Real q = s.dot(output_d) / s.dot(s);
                output_d = output_d - q * s;
                Real mean_error = output_d.array().abs().mean();
                Vector p = (input_d * s) / (s.dot(s));
                // input_d' = input_d - p * s % minimal
                input_d = input_d - p * s.transpose();
 
#ifdef SEPARATE_REGRESSION_RECONSTRUCTION
                projections[d].push_back(Projection(p, u, q));
#endif
                // Step 2: Parameter reconstruction (uses u,p,q from first step)
                // For separation these projection parameters can be stored in a std::vector
                s = identity.transpose() * u;
                parameters.row(d) += q * s.transpose();
                identity -= p * s.transpose();
 
                // Step 3: Stopping condition
                if (mean_error < threshold)
                {
                    break;
                }
            } // for (int k)
        } // for (int d)
 
        return parameters;
    }
 
 
}