卤化物中的 Cholesky 分解答案

【问题标题】：Cholesky decomposition in Halide卤化物中的 Cholesky 分解
【发布时间】：2023-03-28 23:00:01
【问题描述】：

我正在尝试在 Halide 中实现 Cholesky 分解。部分常见算法（如 crout）包括对三角矩阵的迭代。在某种程度上，分解的对角元素是通过从输入矩阵的对角元素中减去部分列和来计算的。列和是根据输入矩阵的三角形部分的平方元素计算的，不包括对角元素。

使用 BLAS，C++ 中的代码如下所示：

double* a; /* input matrix */
int n; /* dimension */
const int c__1 = 1;
const double c_b12 = 1.;
const double c_b10 = -1.;

for (int j = 0; j < n; ++j) {
    double ajj = a[j + j * n] - ddot(&j, &a[j + n], &n, &a[j + n], &n);
    ajj = sqrt(ajj);
    a[j + j * n] = ajj;
    if (j < n) {
            int i__2 = n - j;
            dgemv("No transpose", &i__2, &j, &c_b10, &a[j + 1 + n], &n, &a[j + n], &b, &c_b12, &a[j + 1 + j * n], &c__1);
            double d__1 = 1. / ajj;
            dscal(&i__2, &d__1, &a[j + 1 + j * n], &c__1);            
    }
}

我的问题是，Halide 是否可以普遍表达这样的模式？如果是这样，它会是什么样子？

【问题讨论】：

你最终是否在 Halide 中实现了这个？

标签： c++ halide

【解决方案1】：

我认为 Andrew 可能有更完整的答案，但为了及时响应，您可以使用 RDom 谓词（通过 RDom::where 引入）来枚举三角形区域（或将它们推广到更多维度）。模式的草图是：

Halide::RDom triangular(0, extent, 0, extent);
triangular.where(triangular.x < triangular.y);

然后在缩减中使用triangular。

【讨论】：

【解决方案2】：

我曾经有一个用 Halide 写的快速 Cholesky。不幸的是我找不到代码。我将外部循环放在 C 中，并编写了一个很好的块面板更新例程，该例程一次在类似于 32 宽的面板上运行。这是在 Halide 进行三角迭代之前，所以也许你现在可以做得更好。

【讨论】：