Многомерные нормальные плотности с помощью матрицы Вудбери инверсии Лемма


Я пытаюсь ускорить многомерной нормальной плотности распределения, воспользовавшись тем, что матрица ковариации имеет вид а + у С У'. Обратная такая матрица может быть вычислена с помощью Вудбери матрицы Лемма , что позволяет нам принимать обратных либо диагонали или гораздо меньших матриц.

Аргументы функции динамически размера Eigen::Matrix объекты. Матрицы A и C имеют диагональ и C значительно меньше, чем А.

typedef Eigen::Matrix< double, Eigen::Dynamic, 1              > Vec;
typedef Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > Mat;
const double inv_sqrt_2pi(0.3989422804014327);
const double log_two_pi (1.83787706640935);

double evalMultivNormWBDA(const Vec &x, const Vec &meanVec, const Vec &A, const Mat &U, const Mat &C, bool log)
{
    Mat Ainv = A.asDiagonal().inverse();
    Mat Cinv = C.inverse();
    Mat I =  Cinv + U.transpose()*Ainv*U;
    Mat SigInv = Ainv - Ainv * U * I.ldlt().solve(U.transpose() * Ainv);
    Eigen::LLT<Mat> lltSigInv(SigInv);
    Mat L = lltSigInv.matrixL(); // LL' = Sig^{-1}
    double quadform = (L * (x-meanVec)).squaredNorm();    
    unsigned int dim = x.rows();
    if (log){

        // calculate log-determinant using cholesky decomposition (assumes symmetric and positive definite)
        double halfld (0.0);
        // add up log of diagnols of Cholesky L
        for(size_t i = 0; i < dim; ++i){
            halfld += std::log(L(i,i));
        }

        return -.5*log_two_pi * dim + halfld - .5*quadform;


    }else{  // not the log density
        double normConst = std::pow(inv_sqrt_2pi, dim) * L.determinant();
        return normConst * std::exp(-.5* quadform);
    }
}

Для 9-мерных векторов это на самом деле медленнее, чем функция у меня, что не воспользоваться леммой. Как я могу сделать это быстрее?

Вот это общее многомерной нормальной функции плотности наряду с некоторыми код грубые сроки. Даже для векторов размерности 100, это медленнее.

double evalMultivNorm(const Vec &x, const Vec &meanVec, const Mat &covMat, bool log)
{
    // from Eigen: Remember that Cholesky decompositions are not rank-revealing. 
    /// This LLT decomposition is only stable on positive definite matrices, 
    // use LDLT instead for the semidefinite case. Also, do not use a Cholesky 
    // decomposition to determine whether a system of equations has a solution.
    Eigen::LLT<Mat> lltM(covMat);
    double quadform = (lltM.matrixL().solve(x-meanVec)).squaredNorm();
    size_t dim = x.rows();
    if (log){

        // calculate log-determinant using cholesky decomposition too
        double ld (0.0);
        Mat L = lltM.matrixL(); // the lower diagonal L such that M = LL^T

        // add up log of diagnols of Cholesky L
        for(size_t i = 0; i < dim; ++i){
            ld += std::log(L(i,i));
        }
        ld *= 2; // covMat = LL^T

        return -.5*log_two_pi * dim - .5*ld - .5*quadform;


    }else{  // not the log density
        double normConst = std::pow(inv_sqrt_2pi, dim) / lltM.matrixL().determinant();
        return normConst * std::exp(-.5* quadform);
    }

}

typedef std::chrono::high_resolution_clock Clock;


int main(int argc, char **argv)
{
    auto begin = Clock::now();


    unsigned int size(100);
    unsigned int numiters(1e5);
    Vec x = Vec::Random(size);
    Vec mean = Vec::Zero(size);
    Mat U = Mat::Random(size,1);
    Mat C = Mat::Identity(1,1);
    Mat Amat = Mat::Identity(size,size) * .01;
    Vec Avec = Vec::Constant(size,.01);
    Mat cov = U * C * U.transpose() + Amat;
    for(size_t i = 0; i < numiters; ++i){
        evalMultivNorm(x, mean, cov, true); // 134 milliseconds with ten, 1588 with 50, 7718  100
        //evalMultivNormWBDA(x, mean, Avec, U, C, true); // 368,  2687 with 50, 15583 with 100
    }






    auto end = Clock::now();
    std::cout << "Delta t2-t1: " 
              << std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count()
              << " milliseconds" << std::endl;
    return 0;
}


182
6
задан 26 февраля 2018 в 10:02 Источник Поделиться
Комментарии