Cholesky with floats 128 (Python)

Question

I am trying to use floats 128 to improve the accuracy of very "rough" simulations. However I have a bottleneck with the NumPy's Cholesky decomposition that accepts no more than floats 64.

Question: 1) Is there a way to do Cholesky with np.float128? (I can do with using another library.)

See below some code to reproduce the problem:

import numpy as np

def get_fractional_std_matrix(t: np.ndarray, H: float):
    s = t[np.newaxis, 1:]
    u = t[1:, np.newaxis]
    H2 = 2 * H
    cov_matrix = ((s ** H2) + (u ** H2) - (np.abs(s - u) ** H2)) / 2

    N = len(t)
    std_matrix = np.zeros(shape=(N, N), dtype=t.dtype)
    std_matrix[1:, 1:] = np.linalg.cholesky(cov_matrix.astype(t.dtype))
    return std_matrix

H = 0.1
t = np.linspace(start=0, stop=1, num=100).astype(np.longdouble)
get_fractional_std_matrix(t=t, H=H)  # TypeError: array type float128 is unsupported in linalg

Funnily enough, the addition of .astype(dtype) inside cholesky gives me almost the same error on a 64 bit only machine:

# TypeError: array type float64 is unsupported in linalg

Answer 1

So I did it "by hand" and it seems to work fine with non-complex types including np.float128.

def cholesky(matrix):
    decomposition = np.zeros_like(matrix)
    for i in range(len(matrix)):
        temp_sum = matrix[i:, i] - np.einsum("jk, k -> j", decomposition[i:, :i], decomposition[i, :i])
        decomposition[i, i] = np.sqrt(temp_sum[0])
        decomposition[i+1:, i] = temp_sum[1:] / decomposition[i, i]
    return decomposition