cholesky.hpp

Go to the documentation of this file.

#pragma once
 
#include "NumCpp/Core/Internal/Error.hpp"
#include "NumCpp/Core/Internal/StaticAsserts.hpp"
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/NdArray.hpp"
 
namespace nc::linalg
{
    //============================================================================
    // Method Description:
    template<typename dtype>
    NdArray<double> cholesky(const NdArray<dtype>& inMatrix)
    {
        STATIC_ASSERT_ARITHMETIC(dtype);
 
        const auto shape = inMatrix.shape();
        if (!shape.issquare())
        {
            THROW_RUNTIME_ERROR("Input matrix should be square.");
        }
 
        auto lMatrix = inMatrix.template astype<double>();
 
        for (uint32 row = 0; row < shape.rows; ++row)
        {
            for (uint32 col = row + 1; col < shape.cols; ++col)
            {
                lMatrix(row, col) = 0.;
            }
        }
 
        for (uint32 k = 0; k < shape.cols; ++k)
        {
            const double& a_kk = lMatrix(k, k);
 
            if (a_kk > 0.)
            {
                lMatrix(k, k) = std::sqrt(a_kk);
 
                for (uint32 i = k + 1; i < shape.rows; ++i)
                {
                    lMatrix(i, k) /= lMatrix(k, k);
 
                    for (uint32 j = k + 1; j <= i; ++j)
                    {
                        lMatrix(i, j) -= lMatrix(i, k) * lMatrix(j, k);
                    }
                }
            }
            else
            {
                THROW_RUNTIME_ERROR("Matrix is not positive definite.");
            }
        }
 
        return lMatrix;
    }
} // namespace nc::linalg