include/maths/common/CLinearAlgebraTools.h (676 lines of code) (raw):
/*
* Copyright Elasticsearch B.V. and/or licensed to Elasticsearch B.V. under one
* or more contributor license agreements. Licensed under the Elastic License
* 2.0 and the following additional limitation. Functionality enabled by the
* files subject to the Elastic License 2.0 may only be used in production when
* invoked by an Elasticsearch process with a license key installed that permits
* use of machine learning features. You may not use this file except in
* compliance with the Elastic License 2.0 and the foregoing additional
* limitation.
*/
#ifndef INCLUDED_ml_maths_common_CLinearAlgebraTools_h
#define INCLUDED_ml_maths_common_CLinearAlgebraTools_h
#include <core/CLogger.h>
#include <core/CStringUtils.h>
#include <maths/common/CLinearAlgebra.h>
#include <maths/common/CLinearAlgebraEigen.h>
#include <maths/common/ImportExport.h>
#include <cmath>
#include <cstddef>
#include <limits>
#include <ostream>
#include <vector>
namespace ml {
namespace maths {
namespace common {
namespace linear_algebra_tools_detail {
struct VectorTag;
struct MatrixTag;
struct VectorVectorTag;
struct MatrixMatrixTag;
struct VectorScalarTag;
struct MatrixScalarTag;
struct ScalarVectorTag;
struct ScalarMatrixTag;
template<typename TAG>
struct SSqrt {};
//! Component-wise sqrt for a vector.
template<>
struct SSqrt<VectorTag> {
template<typename VECTOR>
static void calculate(std::size_t d, VECTOR& result) {
for (std::size_t i = 0; i < d; ++i) {
result(i) = std::sqrt(result(i));
}
}
};
//! Element-wise sqrt for a symmetric matrix.
template<>
struct SSqrt<MatrixTag> {
template<typename MATRIX>
static void calculate(std::size_t d, MATRIX& result) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
result(i, j) = std::sqrt(result(i, j));
}
}
}
};
template<typename TAG>
struct SMin {};
//! Component-wise minimum for a vector.
template<>
struct SMin<VectorVectorTag> {
template<typename VECTOR>
static void calculate(std::size_t d, const VECTOR& lhs, VECTOR& rhs) {
for (std::size_t i = 0; i < d; ++i) {
rhs(i) = std::min(lhs(i), rhs(i));
}
}
};
//! Component-wise minimum for a vector.
template<>
struct SMin<VectorScalarTag> {
template<typename VECTOR, typename T>
static void calculate(std::size_t d, VECTOR& lhs, const T& rhs) {
for (std::size_t i = 0; i < d; ++i) {
lhs(i) = std::min(lhs(i), rhs);
}
}
};
//! Component-wise minimum for a vector.
template<>
struct SMin<ScalarVectorTag> {
template<typename T, typename VECTOR>
static void calculate(std::size_t d, const T& lhs, VECTOR& rhs) {
for (std::size_t i = 0; i < d; ++i) {
rhs(i) = std::min(rhs(i), lhs);
}
}
};
//! Element-wise minimum for a symmetric matrix.
template<>
struct SMin<MatrixMatrixTag> {
template<typename MATRIX>
static void calculate(std::size_t d, const MATRIX& lhs, MATRIX& rhs) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
rhs(i, j) = std::min(lhs(i, j), rhs(i, j));
}
}
}
};
//! Element-wise minimum for a symmetric matrix.
template<>
struct SMin<MatrixScalarTag> {
template<typename MATRIX, typename T>
static void calculate(std::size_t d, MATRIX& lhs, const T& rhs) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
lhs(i, j) = std::min(lhs(i, j), rhs);
}
}
}
};
//! Element-wise minimum for a symmetric matrix.
template<>
struct SMin<ScalarMatrixTag> {
template<typename T, typename MATRIX>
static void calculate(std::size_t d, const T& lhs, MATRIX& rhs) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
rhs(i, j) = std::min(lhs, rhs(i, j));
}
}
}
};
template<typename TAG>
struct SMax {};
//! Component-wise maximum for a vector.
template<>
struct SMax<VectorVectorTag> {
template<typename VECTOR>
static void calculate(std::size_t d, const VECTOR& lhs, VECTOR& rhs) {
for (std::size_t i = 0; i < d; ++i) {
rhs(i) = std::max(lhs(i), rhs(i));
}
}
};
//! Component-wise maximum for a vector.
template<>
struct SMax<VectorScalarTag> {
template<typename VECTOR, typename T>
static void calculate(std::size_t d, VECTOR& lhs, const T& rhs) {
for (std::size_t i = 0; i < d; ++i) {
lhs(i) = std::max(lhs(i), rhs);
}
}
};
//! Component-wise maximum for a vector.
template<>
struct SMax<ScalarVectorTag> {
template<typename T, typename VECTOR>
static void calculate(std::size_t d, const T& lhs, VECTOR& rhs) {
for (std::size_t i = 0; i < d; ++i) {
rhs(i) = std::max(rhs(i), lhs);
}
}
};
//! Element-wise maximum for a symmetric matrix.
template<>
struct SMax<MatrixMatrixTag> {
template<typename MATRIX>
static void calculate(std::size_t d, const MATRIX& lhs, MATRIX& rhs) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
rhs(i, j) = std::max(lhs(i, j), rhs(i, j));
}
}
}
};
//! Element-wise maximum for a symmetric matrix.
template<>
struct SMax<MatrixScalarTag> {
template<typename MATRIX, typename T>
static void calculate(std::size_t d, MATRIX& lhs, const T& rhs) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
lhs(i, j) = std::max(lhs(i, j), rhs);
}
}
}
};
//! Element-wise maximum for a symmetric matrix.
template<>
struct SMax<ScalarMatrixTag> {
template<typename T, typename MATRIX>
static void calculate(std::size_t d, const T& lhs, MATRIX& rhs) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
rhs(i, j) = std::max(lhs, rhs(i, j));
}
}
}
};
template<typename TAG>
struct SFabs {};
//! Component-wise fabs for a vector.
template<>
struct SFabs<VectorTag> {
template<typename VECTOR>
static void calculate(std::size_t d, VECTOR& result) {
for (std::size_t i = 0; i < d; ++i) {
result(i) = std::fabs(result(i));
}
}
};
//! Element-wise fabs for a symmetric matrix.
template<>
struct SFabs<MatrixTag> {
template<typename MATRIX>
static void calculate(std::size_t d, MATRIX& result) {
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j <= i; ++j) {
result(i, j) = std::fabs(result(i, j));
}
}
}
};
#define INVERSE_QUADRATIC_PRODUCT(T, N) \
MATHS_COMMON_EXPORT \
maths_t::EFloatingPointErrorStatus inverseQuadraticProduct( \
std::size_t d, const CSymmetricMatrixNxN<T, N>& covariance, \
const CVectorNx1<T, N>& residual, double& result, bool ignoreSingularSubspace)
INVERSE_QUADRATIC_PRODUCT(CFloatStorage, 2);
INVERSE_QUADRATIC_PRODUCT(CFloatStorage, 3);
INVERSE_QUADRATIC_PRODUCT(CFloatStorage, 4);
INVERSE_QUADRATIC_PRODUCT(CFloatStorage, 5);
INVERSE_QUADRATIC_PRODUCT(double, 2);
INVERSE_QUADRATIC_PRODUCT(double, 3);
INVERSE_QUADRATIC_PRODUCT(double, 4);
INVERSE_QUADRATIC_PRODUCT(double, 5);
#undef INVERSE_QUADRATIC_PRODUCT
MATHS_COMMON_EXPORT
maths_t::EFloatingPointErrorStatus
inverseQuadraticProduct(std::size_t d,
const CSymmetricMatrix<CFloatStorage>& covariance,
const CVector<CFloatStorage>& residual,
double& result,
bool ignoreSingularSubspace);
MATHS_COMMON_EXPORT
maths_t::EFloatingPointErrorStatus
inverseQuadraticProduct(std::size_t d,
const CSymmetricMatrix<double>& covariance,
const CVector<double>& residual,
double& result,
bool ignoreSingularSubspace);
#define GAUSSIAN_LOG_LIKELIHOOD(T, N) \
MATHS_COMMON_EXPORT \
maths_t::EFloatingPointErrorStatus gaussianLogLikelihood( \
std::size_t d, const CSymmetricMatrixNxN<T, N>& covariance, \
const CVectorNx1<T, N>& residual, double& result, bool ignoreSingularSubspace)
GAUSSIAN_LOG_LIKELIHOOD(CFloatStorage, 2);
GAUSSIAN_LOG_LIKELIHOOD(CFloatStorage, 3);
GAUSSIAN_LOG_LIKELIHOOD(CFloatStorage, 4);
GAUSSIAN_LOG_LIKELIHOOD(CFloatStorage, 5);
GAUSSIAN_LOG_LIKELIHOOD(double, 2);
GAUSSIAN_LOG_LIKELIHOOD(double, 3);
GAUSSIAN_LOG_LIKELIHOOD(double, 4);
GAUSSIAN_LOG_LIKELIHOOD(double, 5);
#undef GAUSSIAN_LOG_LIKELIHOOD
MATHS_COMMON_EXPORT
maths_t::EFloatingPointErrorStatus
gaussianLogLikelihood(std::size_t d,
const CSymmetricMatrix<CFloatStorage>& covariance,
const CVector<CFloatStorage>& residual,
double& result,
bool ignoreSingularSubspace);
MATHS_COMMON_EXPORT
maths_t::EFloatingPointErrorStatus
gaussianLogLikelihood(std::size_t d,
const CSymmetricMatrix<double>& covariance,
const CVector<double>& residual,
double& result,
bool ignoreSingularSubspace);
//! Shared implementation of Gaussian sampling.
#define SAMPLE_GAUSSIAN(T, N) \
MATHS_COMMON_EXPORT \
void sampleGaussian(std::size_t n, const CVectorNx1<T, N>& mean, \
const CSymmetricMatrixNxN<T, N>& covariance, \
std::vector<CVectorNx1<double, N>>& result)
SAMPLE_GAUSSIAN(CFloatStorage, 2);
SAMPLE_GAUSSIAN(CFloatStorage, 3);
SAMPLE_GAUSSIAN(CFloatStorage, 4);
SAMPLE_GAUSSIAN(CFloatStorage, 5);
SAMPLE_GAUSSIAN(double, 2);
SAMPLE_GAUSSIAN(double, 3);
SAMPLE_GAUSSIAN(double, 4);
SAMPLE_GAUSSIAN(double, 5);
#undef SAMPLE_GAUSSIAN
MATHS_COMMON_EXPORT
void sampleGaussian(std::size_t n,
const CVector<CFloatStorage>& mean,
const CSymmetricMatrix<CFloatStorage>& covariance,
std::vector<CVector<double>>& result);
MATHS_COMMON_EXPORT
void sampleGaussian(std::size_t n,
const CVector<double>& mean,
const CSymmetricMatrix<double>& covariance,
std::vector<CVector<double>>& result);
//! Shared implementation of the log-determinant function.
#define LOG_DETERMINANT(T, N) \
MATHS_COMMON_EXPORT \
maths_t::EFloatingPointErrorStatus logDeterminant( \
std::size_t d, const CSymmetricMatrixNxN<T, N>& matrix, \
double& result, bool ignoreSingularSubspace)
LOG_DETERMINANT(CFloatStorage, 2);
LOG_DETERMINANT(CFloatStorage, 3);
LOG_DETERMINANT(CFloatStorage, 4);
LOG_DETERMINANT(CFloatStorage, 5);
LOG_DETERMINANT(double, 2);
LOG_DETERMINANT(double, 3);
LOG_DETERMINANT(double, 4);
LOG_DETERMINANT(double, 5);
#undef LOG_DETERMINANT
MATHS_COMMON_EXPORT
maths_t::EFloatingPointErrorStatus logDeterminant(std::size_t d,
const CSymmetricMatrix<CFloatStorage>& matrix,
double& result,
bool ignoreSingularSubspace);
MATHS_COMMON_EXPORT
maths_t::EFloatingPointErrorStatus logDeterminant(std::size_t d,
const CSymmetricMatrix<double>& matrix,
double& result,
bool ignoreSingularSubspace);
}
//! Output for debug.
template<typename T>
std::ostream& operator<<(std::ostream& o, const CSymmetricMatrix<T>& m) {
for (std::size_t i = 0; i < m.rows(); ++i) {
o << "\n ";
for (std::size_t j = 0; j < m.columns(); ++j) {
std::string element = core::CStringUtils::typeToStringPretty(m(i, j));
o << element << std::string(15 - element.size(), ' ');
}
}
return o;
}
//! Output for debug.
template<typename T, std::size_t N>
std::ostream& operator<<(std::ostream& o, const CSymmetricMatrixNxN<T, N>& m) {
for (std::size_t i = 0; i < N; ++i) {
o << "\n ";
for (std::size_t j = 0; j < N; ++j) {
std::string element = core::CStringUtils::typeToStringPretty(m(i, j));
o << element << std::string(15 - element.size(), ' ');
}
}
return o;
}
//! Output for debug.
template<typename T, std::size_t N>
std::ostream& operator<<(std::ostream& o, const CVectorNx1<T, N>& v) {
return o << core::CContainerPrinter::print(v.begin(), v.end());
}
//! Output for debug.
template<typename T>
std::ostream& operator<<(std::ostream& o, const CVector<T>& v) {
return o << core::CContainerPrinter::print(v.begin(), v.end());
}
//! Overload sqrt for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> sqrt(const CVectorNx1<T, N>& v) {
CVectorNx1<T, N> result(v);
linear_algebra_tools_detail::SSqrt<linear_algebra_tools_detail::VectorTag>::calculate(
N, result);
return result;
}
//! Overload sqrt for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N> sqrt(const CSymmetricMatrixNxN<T, N>& m) {
CSymmetricMatrixNxN<T, N> result(m);
linear_algebra_tools_detail::SSqrt<linear_algebra_tools_detail::MatrixTag>::calculate(
N, result);
return result;
}
//! Overload minimum for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> min(const CVectorNx1<T, N>& lhs, const CVectorNx1<T, N>& rhs) {
CVectorNx1<T, N> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::VectorVectorTag>::calculate(
N, lhs, result);
return result;
}
//! Overload minimum for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> min(const CVectorNx1<T, N>& lhs, const T& rhs) {
CVectorNx1<T, N> result(lhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::VectorScalarTag>::calculate(
N, result, rhs);
return result;
}
//! Overload minimum for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> min(const T& lhs, const CVectorNx1<T, N>& rhs) {
CVectorNx1<T, N> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::ScalarVectorTag>::calculate(
N, lhs, result);
return result;
}
//! Overload minimum for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N>
min(const CSymmetricMatrixNxN<T, N>& lhs, const CSymmetricMatrixNxN<T, N>& rhs) {
CSymmetricMatrixNxN<T, N> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::MatrixMatrixTag>::calculate(
N, lhs, result);
return result;
}
//! Overload minimum for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N> min(const CSymmetricMatrixNxN<T, N>& lhs, const T& rhs) {
CSymmetricMatrixNxN<T, N> result(lhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::MatrixScalarTag>::calculate(
N, result, rhs);
return result;
}
//! Overload minimum for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N> min(const T& lhs, const CSymmetricMatrixNxN<T, N>& rhs) {
CSymmetricMatrixNxN<T, N> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::ScalarMatrixTag>::calculate(
N, lhs, result);
return result;
}
//! Overload maximum for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> max(const CVectorNx1<T, N>& lhs, const CVectorNx1<T, N>& rhs) {
CVectorNx1<T, N> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::VectorVectorTag>::calculate(
N, lhs, result);
return result;
}
//! Overload maximum for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> max(const CVectorNx1<T, N>& lhs, const T& rhs) {
CVectorNx1<T, N> result(lhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::VectorScalarTag>::calculate(
N, result, rhs);
return result;
}
//! Overload maximum for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> max(const T& lhs, const CVectorNx1<T, N>& rhs) {
CVectorNx1<T, N> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::ScalarVectorTag>::calculate(
N, lhs, result);
return result;
}
//! Overload maximum for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N>
max(const CSymmetricMatrixNxN<T, N>& lhs, const CSymmetricMatrixNxN<T, N>& rhs) {
CSymmetricMatrixNxN<T, N> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::MatrixMatrixTag>::calculate(
N, lhs, result);
return result;
}
//! Overload maximum for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N> max(const CSymmetricMatrixNxN<T, N>& lhs, const T& rhs) {
CSymmetricMatrixNxN<T, N> result(lhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::MatrixScalarTag>::calculate(
N, result, rhs);
return result;
}
//! Overload maximum for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N> max(const T& lhs, const CSymmetricMatrixNxN<T, N>& rhs) {
CSymmetricMatrixNxN<T, N> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::ScalarMatrixTag>::calculate(
N, lhs, result);
return result;
}
//! Overload std::fabs for CVectorNx1.
template<typename T, std::size_t N>
CVectorNx1<T, N> fabs(const CVectorNx1<T, N>& v) {
CVectorNx1<T, N> result(v);
linear_algebra_tools_detail::SFabs<linear_algebra_tools_detail::VectorTag>::calculate(
N, result);
return result;
}
//! Overload std::fabs for CSymmetricMatrixNxN.
template<typename T, std::size_t N>
CSymmetricMatrixNxN<T, N> fabs(const CSymmetricMatrixNxN<T, N>& m) {
CSymmetricMatrixNxN<T, N> result(m);
linear_algebra_tools_detail::SFabs<linear_algebra_tools_detail::MatrixTag>::calculate(
N, result);
return result;
}
//! Overload sqrt for CVector.
template<typename T>
CVector<T> sqrt(const CVector<T>& v) {
CVector<T> result(v);
linear_algebra_tools_detail::SSqrt<linear_algebra_tools_detail::VectorTag>::calculate(
result.dimension(), result);
return result;
}
//! Overload sqrt for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> sqrt(const CSymmetricMatrix<T>& m) {
CSymmetricMatrix<T> result(m);
linear_algebra_tools_detail::SSqrt<linear_algebra_tools_detail::MatrixTag>::calculate(
result.rows(), result);
return result;
}
//! Overload minimum for CVector.
template<typename T>
CVector<T> min(const CVector<T>& lhs, const CVector<T>& rhs) {
CVector<T> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::VectorVectorTag>::calculate(
result.dimension(), lhs, result);
return result;
}
//! Overload minimum for CVector.
template<typename T>
CVector<T> min(const CVector<T>& lhs, const T& rhs) {
CVector<T> result(lhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::VectorScalarTag>::calculate(
result.dimension(), result, rhs);
return result;
}
//! Overload minimum for CVector.
template<typename T>
CVector<T> min(const T& lhs, const CVector<T>& rhs) {
CVector<T> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::ScalarVectorTag>::calculate(
result.dimension(), lhs, result);
return result;
}
//! Overload minimum for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> min(const CSymmetricMatrix<T>& lhs, const CSymmetricMatrix<T>& rhs) {
CSymmetricMatrix<T> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::MatrixMatrixTag>::calculate(
result.rows(), lhs, result);
return result;
}
//! Overload minimum for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> min(const CSymmetricMatrix<T>& lhs, const T& rhs) {
CSymmetricMatrix<T> result(lhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::MatrixScalarTag>::calculate(
result.rows(), result, rhs);
return result;
}
//! Overload minimum for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> min(const T& lhs, const CSymmetricMatrix<T>& rhs) {
CSymmetricMatrix<T> result(rhs);
linear_algebra_tools_detail::SMin<linear_algebra_tools_detail::ScalarMatrixTag>::calculate(
result.rows(), lhs, result);
return result;
}
//! Overload maximum for CVector.
template<typename T>
CVector<T> max(const CVector<T>& lhs, const CVector<T>& rhs) {
CVector<T> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::VectorVectorTag>::calculate(
result.dimension(), lhs, result);
return result;
}
//! Overload maximum for CVector.
template<typename T>
CVector<T> max(const CVector<T>& lhs, const T& rhs) {
CVector<T> result(lhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::VectorScalarTag>::calculate(
result.dimension(), result, rhs);
return result;
}
//! Overload maximum for CVector.
template<typename T>
CVector<T> max(const T& lhs, const CVector<T>& rhs) {
CVector<T> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::ScalarVectorTag>::calculate(
result.dimension(), lhs, result);
return result;
}
//! Overload maximum for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> max(const CSymmetricMatrix<T>& lhs, const CSymmetricMatrix<T>& rhs) {
CSymmetricMatrix<T> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::MatrixMatrixTag>::calculate(
result.rows(), lhs, result);
return result;
}
//! Overload maximum for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> max(const CSymmetricMatrix<T>& lhs, const T& rhs) {
CSymmetricMatrix<T> result(lhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::MatrixScalarTag>::calculate(
result.rows(), result, rhs);
return result;
}
//! Overload maximum for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> max(const T& lhs, const CSymmetricMatrix<T>& rhs) {
CSymmetricMatrix<T> result(rhs);
linear_algebra_tools_detail::SMax<linear_algebra_tools_detail::ScalarMatrixTag>::calculate(
result.rows(), lhs, result);
return result;
}
//! Overload std::fabs for CVector.
template<typename T>
CVector<T> fabs(const CVector<T>& v) {
CVector<T> result(v);
linear_algebra_tools_detail::SFabs<linear_algebra_tools_detail::VectorTag>::calculate(
result.dimension(), result);
return result;
}
//! Overload std::fabs for CSymmetricMatrix.
template<typename T>
CSymmetricMatrix<T> fabs(const CSymmetricMatrix<T>& m) {
CSymmetricMatrix<T> result(m);
linear_algebra_tools_detail::SFabs<linear_algebra_tools_detail::MatrixTag>::calculate(
result.dimension(), result);
return result;
}
//! Scale the \p i'th row and column by \p scale.
template<typename T, typename MATRIX>
void scaleCovariances(std::size_t i, T scale, MATRIX& m) {
scale = std::sqrt(scale);
for (std::size_t j = 0; j < m.columns(); ++j) {
if (i != j) {
m(i, j) *= scale;
} else {
m(i, i) *= scale * scale;
}
}
}
//! Scale the rows and columns by \p scale.
template<typename VECTOR, typename MATRIX>
void scaleCovariances(const VECTOR& scale, MATRIX& m) {
for (std::size_t i = 0; i < scale.dimension(); ++i) {
scaleCovariances(i, scale(i), m);
}
}
//! Scale the \p i'th row and column by \p scale.
template<typename SCALAR>
void scaleCovariances(std::ptrdiff_t i, SCALAR scale, CDenseMatrix<SCALAR>& m) {
scale = std::sqrt(scale);
for (std::ptrdiff_t j = 0; j < m.cols(); ++j) {
if (i != j) {
m(i, j) = m(j, i) = scale * m(i, j);
} else {
m(i, i) *= scale * scale;
}
}
}
//! Scale the rows and columns by \p scale.
template<typename SCALAR>
void scaleCovariances(const CDenseVector<SCALAR>& scale, CDenseMatrix<SCALAR>& m) {
for (std::ptrdiff_t i = 0; i < scale.size(); ++i) {
scaleCovariances(i, scale(i), m);
}
}
//! Compute the inverse quadratic form \f$x^tC^{-1}x\f$.
//!
//! \param[in] covariance The matrix.
//! \param[in] residual The vector.
//! \param[out] result Filled in with the log likelihood.
//! \param[in] ignoreSingularSubspace If true then we ignore the
//! residual on a singular subspace of m. Otherwise the result is
//! minus infinity in this case.
template<typename T, std::size_t N>
maths_t::EFloatingPointErrorStatus
inverseQuadraticForm(const CSymmetricMatrixNxN<T, N>& covariance,
const CVectorNx1<T, N>& residual,
double& result,
bool ignoreSingularSubspace = true) {
return linear_algebra_tools_detail::inverseQuadraticProduct(
N, covariance, residual, result, ignoreSingularSubspace);
}
//! Compute the log-likelihood for the residual \p x and covariance
//! matrix \p m.
//!
//! \param[in] covariance The matrix.
//! \param[in] residual The vector.
//! \param[out] result Filled in with the log likelihood.
//! \param[in] ignoreSingularSubspace If true then we ignore the
//! residual on a singular subspace of m. Otherwise the result is
//! minus infinity in this case.
template<typename T, std::size_t N>
maths_t::EFloatingPointErrorStatus
gaussianLogLikelihood(const CSymmetricMatrixNxN<T, N>& covariance,
const CVectorNx1<T, N>& residual,
double& result,
bool ignoreSingularSubspace = true) {
return linear_algebra_tools_detail::gaussianLogLikelihood(
N, covariance, residual, result, ignoreSingularSubspace);
}
//! Sample from a Gaussian with \p mean and \p covariance in such
//! a way as to preserve the mean, covariance matrix and some of
//! the quantiles of the generalised c.d.f.
//!
//! \param[in] n The desired number of samples.
//! \param[in] mean The mean of the Gaussian.
//! \param[in] covariance The covariance matrix of the Gaussian.
//! \param[out] result Filled in with the samples.
template<typename T, typename U, std::size_t N>
void sampleGaussian(std::size_t n,
const CVectorNx1<T, N>& mean,
const CSymmetricMatrixNxN<T, N>& covariance,
std::vector<CVectorNx1<U, N>>& result) {
return linear_algebra_tools_detail::sampleGaussian(n, mean, covariance, result);
}
//! Compute the log-determinant of the symmetric matrix \p m.
//!
//! \param[in] matrix The matrix.
//! \param[in] ignoreSingularSubspace If true then we ignore any
//! singular subspace of m. Otherwise, the result is minus infinity.
template<typename T, std::size_t N>
maths_t::EFloatingPointErrorStatus logDeterminant(const CSymmetricMatrixNxN<T, N>& matrix,
double& result,
bool ignoreSingularSubspace = true) {
return linear_algebra_tools_detail::logDeterminant(N, matrix, result, ignoreSingularSubspace);
}
//! Compute the inverse quadratic form \f$x^tC^{-1}x\f$.
//!
//! \param[in] covariance The matrix.
//! \param[in] residual The vector.
//! \param[out] result Filled in with the log likelihood.
//! \param[in] ignoreSingularSubspace If true then we ignore the
//! residual on a singular subspace of m. Otherwise the result is
//! minus infinity in this case.
template<typename T>
maths_t::EFloatingPointErrorStatus
inverseQuadraticForm(const CSymmetricMatrix<T>& covariance,
const CVector<T>& residual,
double& result,
bool ignoreSingularSubspace = true) {
return linear_algebra_tools_detail::inverseQuadraticProduct(
covariance.rows(), covariance, residual, result, ignoreSingularSubspace);
}
//! Compute the log-likelihood for the residual \p x and covariance
//! matrix \p m.
//!
//! \param[in] covariance The covariance matrix.
//! \param[in] residual The residual, i.e. x - mean.
//! \param[out] result Filled in with the log likelihood.
//! \param[in] ignoreSingularSubspace If true then we ignore the
//! residual on a singular subspace of m. Otherwise the result is
//! minus infinity in this case.
template<typename T>
maths_t::EFloatingPointErrorStatus
gaussianLogLikelihood(const CSymmetricMatrix<T>& covariance,
const CVector<T>& residual,
double& result,
bool ignoreSingularSubspace = true) {
return linear_algebra_tools_detail::gaussianLogLikelihood(
covariance.rows(), covariance, residual, result, ignoreSingularSubspace);
}
//! Sample from a Gaussian with \p mean and \p covariance in such
//! a way as to preserve the mean, covariance matrix and some of
//! the quantiles of the generalised c.d.f.
//!
//! \param[in] n The desired number of samples.
//! \param[in] mean The mean of the Gaussian.
//! \param[in] covariance The covariance matrix of the Gaussian.
//! \param[out] result Filled in with the samples.
template<typename T, typename U>
void sampleGaussian(std::size_t n,
const CVector<T>& mean,
const CSymmetricMatrix<T>& covariance,
std::vector<CVector<U>>& result) {
return linear_algebra_tools_detail::sampleGaussian(n, mean, covariance, result);
}
//! Compute the log-determinant of the symmetric matrix \p m.
//!
//! \param[in] matrix The matrix.
//! \param[in] ignoreSingularSubspace If true then we ignore any
//! singular subspace of m. Otherwise, the result is minus infinity.
template<typename T>
maths_t::EFloatingPointErrorStatus logDeterminant(const CSymmetricMatrix<T>& matrix,
double& result,
bool ignoreSingularSubspace = true) {
return linear_algebra_tools_detail::logDeterminant(matrix.rows(), matrix, result,
ignoreSingularSubspace);
}
//! Project the matrix on to \p subspace.
template<typename MATRIX>
inline Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
projectedMatrix(const std::vector<std::size_t>& subspace, const MATRIX& matrix) {
std::size_t d = subspace.size();
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> result(d, d);
for (std::size_t i = 0; i < d; ++i) {
for (std::size_t j = 0; j < d; ++j) {
result(i, j) = matrix(subspace[i], subspace[j]);
}
}
return result;
}
//! Project the vector on to \p subspace.
template<typename VECTOR>
inline Eigen::Matrix<double, Eigen::Dynamic, 1>
projectedVector(const std::vector<std::size_t>& subspace, const VECTOR& vector) {
std::size_t d = subspace.size();
Eigen::Matrix<double, Eigen::Dynamic, 1> result(d);
for (std::size_t i = 0; i < d; ++i) {
result(i) = vector(subspace[i]);
}
return result;
}
}
}
}
#endif // INCLUDED_ml_maths_common_CLinearAlgebraTools_h