GGML PCA

主成分分析

主成分分析（principal components analysis）是一类非常常见的算法，广泛应用于各种需要信息压缩的场景。虽然《深度学习》一书称之为“简单的机器学习算法”，但是PCA通常可以通过确定的数值算法实现。概略来说，PCA是找到一个映射算法，将高维空间的样本集合在尽量不损失信息的前提下，映射到更低的维度，从N维空间映射到M维空间，可以表示为一个 $M*N$的矩阵。PCA的目的就是通过分析样本集合，解出这个映射矩阵。这个求解过程主要有这样几个步骤：

将样本集合表达为$样本数特征数$（$AF$）的矩阵
将这个矩阵规范化，使后续步骤中我们可以得到一个对角线为1的协方差矩阵，这需要先得到无偏的标准差：
- 我们先求得样本集的均值向量，均值向量的每个维度都是所有样本在该特征维度的值的平均
- 用样本集减去这个均值向量，得到一组中心化的数据$B$
- 对这个集合求平方，然后求和，再除以$A-1$得到无偏方差
- 开方得到无偏标准差(Unbiased standard deviation)即样本标准差(Samples standard deviation) 。
- 为了避免样本标准差的某个维度为0，导致后续计算出错，可以给它加上一个极小的正数作为补偿。
- 用B矩阵除以标准差
得到规范化的矩阵X之后，我们计算 $X^TX$
将得到的矩阵除以$A-1$，得到$F*F$的协方差矩阵
求解协方差矩阵的奇异值矩阵$S$和右奇异值$VT$向量矩阵，按从奇异值大到小，取对应的前M个向量，组成PCA转换矩阵使用这个矩阵，对样本做矩阵乘法时，得到的就是降维后的结果向量。
实现

我最初做PCA算法的动机，来自于将 4096 维度的 ollama 嵌入向量降维到 2000 维以下——我的目标是256，那么在这个计算过程中，如果假设样本集数量为一万，计算过程中就会涉及好几次千万级的浮点数乘法。为了加速计算过程，使之可以优化到足以放到 PostgreSQL 内部使用，我使用AI框架 ggml 实现了一个 PCA 方法。 GGML 是一个 C 库，它支持 CUDA、MPX 等硬件加速方案，可以快速处理大规模的矩阵算法。 GGML 的tensor类型支持最多4个张量维度，PCA仅涉及矩阵和向量计算，只需要用到最多两位维度。所以这里仅使用GGML TENSOR的前两个维度。在前述步骤中，基础的矩阵算法如乘法，加法等都有对应的 GGML 算子，但是求解 SVD时，要涉及对矩阵方程的数值解运算，GGML中还没有对应的实现，这里我暂时采用了 lapack 的 SVD 功能。 GGML 本身是惰性计算，先初始化context，然后构造执行图，在计算执行图的时候才进行实际的计算。显然这要求我们在求解 SVD 之前先对执行图求值。虽然执行图本身是可再入的，但是我暂时保守的将pca函数设计为封闭的算式，它不要求调用方先准备好执行图，而是在函数执行过程中自己定义一个执行图。其实这个设计并不完美，GGML没有提供显式释放 graph 的操作，它只能随着ctx的释放而释放。这样一个隐式 graph 每次会造成 2k 左右的泄露——直至我们释放 context 。考虑到每一次计算任务可以使用单独的 context，这个问题也不算大。当然，最理想的是，将来做一个svd算子，这样我们就可以将pca做成和其它 tensor dancer 算子一样，遵循 ggml 的规范。是的，在开发pca算法的过程中，我顺便实现了方差、标准差，无偏方差，无偏标准差，协方差等方法，这些方法的调用结果都是 ggml_tensor * ，这样我可以将它们视作是 GGML 的组合算子。它们都是惰性的，而且会在分解为原子的 GGML 算子后执行。

另外，我也提供了 GGML tensor 的 svd 求解封装——内部调用的lapack，所以它不是惰性的，也无需规划执行图，并且，它假定传入的张量是一个 float 32 矩阵。在这个基础上，pca可以简单的封装为：

    struct ggml_tensor *pca = dancer_pca_force(ctx, samples, down_to, 8);

虽然最后求解svd只能在lapack中进行，无法享受到GGML现代化的GPU加速能力，但是考虑到lapack需要处理的部分仅和样本特征相关，即使样本量很大，这个矩阵参数也不会超过特征数的平方，做为第一版，还是可以接受的。具有了这些统计计算函数

int lapack_svd_rf32(ggml_tensor *tensor_A, ggml_tensor *tensor_S, ggml_tensor *tensor_VT);

// centralization
ggml_tensor *dancer_centralization(ggml_context *ctx, ggml_tensor *tensor_X);

// variance
ggml_tensor *dancer_variance(ggml_context *ctx, ggml_tensor *tensor_X);

// normalize
ggml_tensor *dancer_normalized(ggml_context *ctx, ggml_tensor *tensor_X);

// unbiased variance
ggml_tensor *dancer_unbiased_variance(ggml_context *ctx, ggml_tensor *tensor_X);

// population standard deviation
ggml_tensor *dancer_psdv(ggml_context *ctx, ggml_tensor *tensor_X);
// unbiased standard deviation
ggml_tensor *dancer_usdv(ggml_context *ctx, ggml_tensor *tensor_X);
// covariance
ggml_tensor* dancer_covariance(ggml_context *ctx, ggml_tensor * tensor_X);

之后，tensor dancer 已经初步成为了一个支持现代硬件加速能力的统计计算工具。在 TENSOR DANCER 的 test目录中，我写了一些程序，用来测试和演示这些工具。另外有一些数据处理工具和对照程序，用python编写，放在scripts目录。

写在最后

我研究了一下 GGML 内置的规范化操作算子 ggml_norm，它使用的是总体标准差，而协方差和pca计算需要无偏标准差，在 TENSOR DANCER里我将两种都实现了出来。在编写对照程序时，我发现 numpy 的 std 方法默认计算的是总体标准差，如果要计算无偏标准差，需要设定ddof参数，例如

# 计算每个特征的无偏标准差
std = np.std(data, axis=0, ddof=1)

另外，尽管我使用的 AI 坚持 sklearn 的 pca 计算会自动做规范化预处理，但是经过调试，我发现并非如此，如果不对数据集做规范化预处理，sklearn 的 pca 计算并不准确。后面我仍然需要实现 pgvector 、 tensor dancer matrix， ggml tensor之间的转换程序，并且将它们封装为 pg 插件。特别是 tensor dancer 的内部计算大量的依赖 ggml tensor，但是显然作为序列化格式，matrix更紧凑，有它的用武之地。如何让它们之间的转换和数据传递更高效，是个需要思考的问题。