如何使用np.float128进行Cholesky分解？（可借助其他库）

问题背景 #

为提升低精度模拟的准确性，尝试使用np.float128（对应np.longdouble在部分平台的实现），但遇到瓶颈：NumPy原生的np.linalg.cholesky仅支持最高float64类型，无法处理float128数组，报错如下：

TypeError: array type float128 is unsupported in linalg

复现代码：

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

问题解答：是否可以用np.float128做Cholesky分解？ #

可以，但不能依赖NumPy原生的np.linalg.cholesky，需要借助其他库或手动实现，以下是三种可行方案：

方案1：使用SciPy的scipy.linalg.cholesky #

若你的SciPy是基于MKL编译（Linux/macOS常见），scipy.linalg.cholesky支持longdouble（即float128）类型的矩阵分解。修改代码如下：

import numpy as np
from scipy.linalg import cholesky

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)
    # 替换为scipy的cholesky
    std_matrix[1:, 1:] = cholesky(cov_matrix.astype(t.dtype), lower=True)
    return std_matrix

H = 0.1
t = np.linspace(start=0, stop=1, num=100).astype(np.longdouble)
result = get_fractional_std_matrix(t=t, H=H)
print(result.dtype)  # 输出float128（取决于平台）

方案2：使用PyTorch（Linux平台） #

PyTorch在Linux平台支持torch.float128类型，可通过张量转换实现分解：

import numpy as np
import torch

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)
    
    # 转成PyTorch float128张量
    cov_tensor = torch.tensor(cov_matrix, dtype=torch.float128)
    chol_tensor = torch.linalg.cholesky(cov_tensor)
    # 转回NumPy数组
    std_matrix[1:, 1:] = chol_tensor.numpy()
    return std_matrix

H = 0.1
t = np.linspace(start=0, stop=1, num=100).astype(np.longdouble)
result = get_fractional_std_matrix(t=t, H=H)

方案3：手动实现Cholesky分解（适合小矩阵） #

对于小规模矩阵（如示例中的99x99），可以手动实现Cholesky的迭代算法，完全基于float128运算：

import numpy as np

def cholesky_float128(matrix):
    n = matrix.shape[0]
    L = np.zeros_like(matrix)
    for i in range(n):
        for j in range(i+1):
            if i == j:
                # 计算对角线元素
                sum_val = np.sum(L[i, :i] ** 2)
                L[i, i] = np.sqrt(matrix[i, i] - sum_val)
            else:
                # 计算非对角线元素
                sum_val = np.sum(L[i, :j] * L[j, :j])
                L[i, j] = (matrix[i, j] - sum_val) / L[j, j]
    return L

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:] = cholesky_float128(cov_matrix.astype(t.dtype))
    return std_matrix

H = 0.1
t = np.linspace(start=0, stop=1, num=100).astype(np.longdouble)
result = get_fractional_std_matrix(t=t, H=H)

注意：手动实现的效率远低于库函数，仅适合小矩阵场景。

内容的提问来源于stack exchange，提问作者Louis-Amand