概率张量分解与函数配准的统一框架：光滑重参数化实战

概率张量分解和函数型数据配准，在实际工程里经常被分成两个完全不同的任务处理。前者处理的是带约束的高维数组分解，后者处理的是曲线对齐问题。但这两个任务在同一个几何对象下可以统一起来：单纯形乘积空间（simplicial product spaces）。在这个空间上做优化，最大难点不是模型复杂度，而是因子、概率向量或 warping 函数都带有约束，直接用梯度下降时容易越过边界、产生不可解释结果，或者陷入投影步骤带来的额外偏差。

光滑重参数化（smooth reparameterization）解决的就是这个问题：把带约束参数映射到无约束欧氏空间，让普通优化器和普通采样器可以直接使用。本文会先解释单纯形乘积空间为什么是概率张量分解和函数配准的共同基础，然后给出两个基于 PyTorch 的最小可运行案例：一个是带概率约束的 CP 张量分解，另一个是用单调 warping 函数完成曲线配准。通过这两个案例，可以看清重参数化变换的选择如何影响收敛、梯度和结果合法性。

1. 单纯形乘积空间：两个看似不同的问题为何共享同一几何结构

1.1 单纯形不是简单约束，而是“概率向量”所在流形

先说最基础的概念。一个 D 维单纯形可以写成：

Delta^{D-1} = { p in R^D : p_i >= 0, sum(p_i) = 1 }

这个集合的每一个元素，本质上就是一个合法的离散概率分布。比如一个三分类问题的概率向量[0.2, 0.5, 0.3]，就落在二维单纯形上。之所以叫“单纯形”，是因为D=3时它在三维空间中是一个三角形平面；D=4时是一个三角锥的类比体。

在很多机器学习模型里，我们不希望模型输出的向量是任意实数，而是希望它满足“非负且和为 1”。这时候，模型输出层的设计、损失函数的选择、优化器更新后的投影方式，都会受到单纯形结构的限制。很多工程问题，比如主题模型、混合模型、分解模型中的归一化因子，本质都是在单纯形上做推理和优化。

容易误解的一点是，单纯形并不是“添加一个归一化层”就完事。归一化只是把结果投影回单纯形，但优化过程中梯度如何穿越边界、投影是否保持收敛方向、概率值接近 0 时数值是否稳定，这些问题都需要额外处理。单纯形是一个低维流形，而不是一个普通的开集或闭集。边界上的点（某些分量等于 0）在概率意义上对应退化分布，在优化意义上对应梯度消失或对数计算失败的高发区。

1.2 乘积空间的意义：多组概率向量联合出现

单纯形乘积空间，是多个单纯形的笛卡尔积：

Delta^{D1-1} x Delta^{D2-1} x ... x Delta^{DM-1}

它描述的是“多个概率向量同时出现”的场景。概率张量分解就是典型例子。

以三阶张量X ∈ R^{I×J×K}为例，CP 分解把它近似成 R 个秩一分量的和：

X_{ijk} ≈ sum_{r=1}^{R} A_{ir} * B_{jr} * C_{kr}

如果进一步要求每个因子矩阵的每一列都是概率向量，也就是说，对于每个 r，A[:,r]、B[:,r]、C[:,r]都落在单纯形上，那么这个张量分解就可以解释为多个离散分布的组合。比如在传感器多通道数据中，每个通道的信号占比可以用一个概率向量表示；在文本主题张量中，每个主题的词分布就是单纯形上的点。

这时候，整个待优化的参数空间就是三个单纯形的乘积：

Delta^{I-1} x Delta^{J-1} x Delta^{K-1}

函数型数据配准也有类似结构，只是表现形式更隐蔽。假设两条曲线f(t)和g(t)是通过同一个过程生成的，但横轴发生错位，比如一个波峰在 0.4 处，另一个在 0.7 处。配准的任务是找到一个横轴变换w(t)，使得f(w(t))与g(t)对齐。如果把曲线在一组离散时间点上的归一化高度看成概率质量，那么配准就等价于在时间轴上移动这些概率质量，使其在两个分布之间建立对应关系。离散化之后，每条曲线在不同区间上的质量比重也是一个概率向量，多个区间拼接起来就是多个单纯形的乘积空间。

因此，概率张量分解中的因子矩阵，和函数配准中的 warping 分布，虽然应用场景不同，但都共享“多个概率向量协同变化”的数学结构。理解了这一点，再看重参数化就顺理成章。

1.3 在这个空间上做优化的直接困难

直接在单纯形乘积空间上做梯度下降，会遇到三类问题。

第一类是边界问题。梯度更新后，参数可能落在单纯形之外，比如某个分量变成负数，或者总和不再是 1。常规做法是更新后做投影或归一化，但投影会破坏梯度的方向，导致目标函数不降反升。

第二类是退化问题。概率向量逼近边界时，比如某个分量趋近 0，基于对数的正则项会迅速变成无穷大。即使不崩溃，梯度也会在边界附近产生不合理的巨大值，让优化过程来回震荡。

第三类是采样和推断问题。在贝叶斯视角下，概率张量分解通常依赖 Dirichlet 先验或 Logistic-Normal 先验。从 Dirichlet 分布采样后再计算梯度，估计出的梯度方差往往较大，影响变分推断的收敛。

重参数化的关键价值，就是避开这些直接在约束空间上处理问题的麻烦。

2. 光滑重参数化：用无约束欧氏空间替代约束空间

2.1 核心思想：定义一枚光滑双射

光滑重参数化的基本策略，是找一个从开单纯形到欧氏空间的光滑双射（同胚），把带约束的变量p映射成无约束变量z，然后在z空间里执行优化或采样。反向变换再把z映射回单纯形。

这样做的好处很直接：

• 不用做投影，因为从z通过光滑映射得到的结果天然满足概率约束；
• 梯度方向更平稳，因为映射本身是可微的；
• 可以复用普通的 Adam、SGD 等优化器，不需要实现投影梯度或 Riemannian 优化器。

这个思路的关键，是映射必须是“光滑”的。如果映射存在不可导点，比如绝对值、阶跃函数，那么反向传播在不可导点附近会得到不稳定的梯度，甚至直接中断。

2.2 softmax 的问题：冗余自由度

最容易想到的映射是 softmax：

p_i = exp(z_i) / sum_j exp(z_j)

softmax 确实能把任意实数向量映射到单纯形内部，而且完全可微。但它有一个工程上容易忽略的问题：softmax 不是单射。对任意常数c，把z_i全部加上c后，softmax(z + c)和softmax(z)完全相等。也就是说，softmax 把无穷多个不同的z映射到同一个p。

这意味着，如果我们直接优化z，模型会存在冗余自由度。对于一个 D 维概率向量，p只有 D-1 个有效自由度，而z有 D 个自由度。冗余本身不一定会让优化失败，但它会带来平坦方向、初始化敏感和参数解释性差的问题。

更规范的做法是对z施加零均值约束，也就是中心化。常用的变换是中心化对数比变换（Centered Log-Ratio，CLR）：

z_i = log(p_i) - (1/D) * sum_j log(p_j)

反变换通过 softmax 完成：

p_i = softmax(z_i)

此时sum(z) = 0，自由度从 D 降到 D-1，映射是开单纯形到零均值子空间的光滑同胚。这个变换在成分数据分析和变分推断中都非常常用。

2.3 为什么“光滑”在变分推断里如此重要

在概率张量分解的贝叶斯推断中，重参数化不只是优化技巧，更是梯度估计的基础。

变分推断通常需要优化一个下界。当变分分布定义在单纯形上时，最直接的选择是 Dirichlet 分布。但 Dirichlet 分布本身在采样和梯度传播上并不友好：直接从 Dirichlet 采样得到样本后，对变分参数求梯度，需要用得分函数估计器，方差通常较大。

一个替代方案是使用 Logistic-Normal 分布。它的做法是：先从一个高斯分布中采样得到z，然后通过p = softmax(z)映射到单纯形。因为高斯采样可以写成“标准正态噪声 + 均值 + 标准差”的形式，梯度可以通过这个确定性路径回传，方差远低于直接估计。这个技巧就是重参数化梯度估计，也就是变分自编码器里常见的“reparameterization trick”。

这个技巧能成立的前提，正是z -> p的映射是光滑的。如果 softmax 换成不可导的映射，梯度就无法通过采样节点回传。

下面是三种常见映射的对比。

2.4 数值稳定性和边界处理

即使使用了光滑映射，数值稳定性仍然需要单独处理。

当p的某个分量非常接近 0 时，log(p)会变成很大的负数，反向传播时梯度也可能爆炸。应对方法通常有三种：

• 在取对数前对p做下限截断，比如torch.clamp(p, min=1e-8)；
• 在反向变换后用一个小阈值防止概率为 0；
• 在损失函数中使用 log-sum-exp 或 log-softmax 形式，避免显式构造中间概率。

学习环境下可以图省事直接截断，但生产环境里要意识到截断阈值会影响梯度尺度。阈值设得太大，概率值被压平，模型的表达能力下降；阈值设得太小，接近边界时仍然可能产生极端梯度。

3. 概率张量分解最小案例：CLR 重参数化 + CP 分解

3.1 问题设定：约束每列都是概率分布

这一节的目标是用 PyTorch 实现一个带概率约束的三阶 CP 张量分解。给定观测张量X ∈ R^{I×J×K}，我们需要找到三个因子矩阵A ∈ R^{I×R}、B ∈ R^{J×R}、C ∈ R^{K×R}，使得：

X_{ijk} ≈ sum_{r=1}^{R} A_{ir} * B_{jr} * C_{kr}

并且对每个r，矩阵的每一列A[:,r]、B[:,r]、C[:,r]都满足“非负且和为 1”。

这个设定在很多场景中是有实际意义的。比如X是一个三维共现计数张量，每个模式表示一个类别变量，那么A[:,r]可以解释为第 r 个隐成分在第一个类别上的条件概率分布。分解结果天然可解释，不需要在分解后手动归一化。

为了在优化过程中保持概率约束，我们不直接优化A、B、C，而是优化它们的 CLR 坐标ZA、ZB、ZC。每一步前向传播时，用ZA -> A、ZB -> B、ZC -> C的映射把无约束坐标转回概率向量。

3.2 环境与依赖

下面代码基于以下环境，版本只需要支持 PyTorch 自动求导即可。

不需要安装额外复杂依赖。学习环境直接使用 CPU 运行，小规模张量训练几十秒即可完成。

3.3 合成数据生成

为了验证分解是否能恢复真实结构，先手工生成一组真实因子矩阵，然后用这些因子构造一个三阶张量并加入噪声。这样我们可以在实验结束后，用重构误差来判断模型是否学到有效信息。

import torch torch.manual_seed(0) I, J, K = 10, 12, 15 R = 4 # 从 Dirichlet 分布生成真实的概率因子 def dirichlet_sample(alpha, size): return torch.distributions.Dirichlet(alpha).sample(size) A_true = dirichlet_sample(torch.ones(I), (R,)).T # shape [I, R] B_true = dirichlet_sample(torch.ones(J), (R,)).T # shape [J, R] C_true = dirichlet_sample(torch.ones(K), (R,)).T # shape [K, R] # 构造无噪声张量 X_true = torch.einsum('ir,jr,kr->ijk', A_true, B_true, C_true) # 加入轻微噪声，模拟观测误差 noise = 0.02 * torch.randn_like(X_true) X = torch.clamp(X_true + noise, min=0.0)

这里使用 Dirichlet 分布生成真实因子，好处是每个真实因子天然满足概率约束。后续优化时，我们用 CLR 坐标去拟合这些概率因子，而不是直接去拟合无约束实数。

3.4 完整实现

下面实现完整的 CLR 变换和训练循环。

import torch import torch.nn.functional as F def clr_from_prob(p, eps=1e-8): p = torch.clamp(p, min=eps) log_p = torch.log(p) return log_p - log_p.mean(dim=0, keepdim=True) def inv_clr(z): return F.softmax(z, dim=0) def cp_reconstruct(A, B, C): return torch.einsum('ir,jr,kr->ijk', A, B, C) def dirichlet_kl(p, alpha_uniform=1.0): p = torch.clamp(p, min=1e-8) k = p.shape[0] uni = torch.full_like(p, 1.0 / k) return (p * (torch.log(p) - torch.log(uni))).sum() # 初始化无约束坐标 ZA = torch.randn(I, R, requires_grad=True) * 0.1 ZB = torch.randn(J, R, requires_grad=True) * 0.1 ZC = torch.randn(K, R, requires_grad=True) * 0.1 optimizer = torch.optim.Adam([ZA, ZB, ZC], lr=0.01) mse_loss_fn = torch.nn.MSELoss() loss_history = [] for step in range(500): optimizer.zero_grad() A = inv_clr(ZA) B = inv_clr(ZB) C = inv_clr(ZC) X_hat = cp_reconstruct(A, B, C) mse_loss = mse_loss_fn(X_hat, X) kl_reg = 0.01 * ( dirichlet_kl(A) + dirichlet_kl(B) + dirichlet_kl(C) ) loss = mse_loss + kl_reg loss.backward() optimizer.step() loss_history.append(loss.item()) if (step + 1) % 100 == 0: print(f"step {step+1:4d}, loss = {loss.item():.6f}, mse = {mse_loss.item():.6f}")

这段代码的关键点可以拆成三个。第一，ZA、ZB、ZC初始化为randn * 0.1，不要初始化为零，否则 CLR 坐标全为 0 时 softmax 会给出均匀分布，梯度更新方向可能不对称。第二，inv_clr使用F.softmax(z, dim=0)，由于 softmax 对整体平移不变，加上 CLR 的中心化约束后，坐标是唯一的。第三，dirichlet_kl加了很小的权重0.01，目的是让因子不要完全退化成边界上的点，对结果影响不大，但能体现概率分解的贝叶斯背景。

3.5 代码关键点逐段解释

为什么不在A空间直接做梯度下降，而要绕一圈经过Z？

直接优化A时，需要保证每列非负且和为 1。即使初始化合法，普通梯度更新一步后就可能破坏约束。常见做法是更新后再做归一化或投影，但投影后的结果并不是梯度下降的最优方向，迭代次数会增加，甚至出现振荡。

而使用Z坐标后，每次前向计算都会生成合法的概率矩阵。优化器只需要在无约束实数空间里移动，不需要关注投影规则。这个模式可以推广到更复杂的模型：只要约束空间能找到一个光滑双射到欧氏空间，就可以用同样的方式做无约束优化。

另一个需要注意的点是因子分解本身的可辨识性。CP 分解中，A、B、C存在列置换和尺度模糊。加上概率约束后，尺度被固定为列和为 1，但仍存在列置换问题。因此验证模型时，更适合看重构张量误差，而不是直接逐列比较A和A_true。

4. 函数型数据配准：warping 函数如何做光滑重参数化

4.1 配准问题定义

函数型数据配准在医学影像、时序信号、运动捕捉、光谱分析里非常常见。比如多个受试者的心电信号，虽然形态相似，但每个波峰出现的时间会有偏移。如果直接进行逐点平均，波峰会因为错位而被抹平。配准的目标是把这些时间偏移对齐。

用数学语言说，给定模板函数f(t)和观测函数g(t)，需要找到一个横轴变换w(t)，使得：

f(w(t)) ≈ g(t)

w(t)需要满足三个基本约束：

• w(t)的值域在时间区间内；
• w(t)严格单调递增；
• 通常还要求w(0)=0、w(1)=1，保持时间边界不变。

如果w(t)不单调，对齐过程会把时间折叠，导致曲线变形出现交叉，这在很多应用中是不可接受的。

4.2 用“非负平方 + 积分归一化”构造单调 warping

直接参数化一个严格单调递增函数并不容易。一个优雅的方法是借助积分。

令g(t)为任意光滑函数，对它取平方：

h(t) = g(t)^2 >= 0

那么h(t)的积分H(t) = ∫_0^t h(s) ds必然单调不减。再归一化：

w(t) = H(t) / H(T)

就可以得到w(0)=0、w(1)=1的单调递增函数。如果h(t)在某个区间上恒为零，w会出现平台，这可以解释为时间在局部被压缩；只要平台不是全部区域，w仍然满足单调不减。

这个构造方式的光滑性取决于基函数的选择。如果g(t)用光滑基函数展开，比如 B 样条或傅里叶基，那么整个w(t)就是光滑的。优化变量变成基函数系数，不再需要担心单调约束。

4.3 最小案例：两条错位高斯峰

下面的例子中，模板函数和观测函数都是由两个高斯峰组成的曲线，观测函数的峰整体向右偏移。目标是学习一个 warping 函数，让模板曲线在 warp 后与观测曲线对齐。

import torch import math torch.manual_seed(1) N = 128 t = torch.linspace(0.0, 1.0, N) def template(t): return ( torch.exp(-((t - 0.35) / 0.08) ** 2) + 0.6 * torch.exp(-((t - 0.70) / 0.10) ** 2) ) def observed(t): return ( torch.exp(-((t - 0.60) / 0.08) ** 2) + 0.6 * torch.exp(-((t - 0.95) / 0.10) ** 2) ) f = template(t) g = observed(t)

观测函数相当于把模板在横轴上向右平移了大约 0.25，然后截断在[0,1]区间内。注意第二个峰从 0.70 移到 0.95，已经非常靠近右边界，因此 warping 函数需要在前半段快速追赶，后半段保持相对平缓。

4.4 优化与损失函数

使用一组傅里叶基函数来参数化g(t)，然后计算 w：

n_basis = 12 basis = [torch.ones(N, 1)] for k in range(1, n_basis): basis.append(torch.cos(2 * math.pi * k * t).reshape(-1, 1)) basis.append(torch.sin(2 * math.pi * k * t).reshape(-1, 1)) B = torch.cat(basis, dim=1) # shape [N, 2*n_basis - 1] coeff = torch.randn(B.shape[1], requires_grad=True) * 0.5 optimizer = torch.optim.Adam([coeff], lr=0.01) for step in range(600): optimizer.zero_grad() u = B @ coeff h = u ** 2 dt = 1.0 / (N - 1) H = torch.cumsum(h, dim=0) * dt w = H / (H[-1] + 1e-8) f_warp = template(w) align_loss = torch.mean((f_warp - g) ** 2) # warping 导数近似，用于平滑正则 w_prime = h / (H[-1] + 1e-8) smooth_loss = torch.mean((w_prime - 1.0) ** 2) lam = 0.1 loss = align_loss + lam * smooth_loss loss.backward() optimizer.step() if (step + 1) % 100 == 0: print(f"step {step+1:4d}, align = {align_loss.item():.6f}, smooth = {smooth_loss.item():.6f}")

这里的关键是，coeff是无约束实数，经过基函数展开、平方、累积积分、归一化后得到的w天然满足单调和端点约束。因此优化器始终在合法函数空间内搜索。

损失函数包含两项。第一项是对齐误差，要求f(w(t))逼近g(t)。第二项是平滑正则，要求w'(t)尽量接近 1，也就是 warping 不要过于剧烈。lam=0.1时，对齐优先，但平滑项能防止过度扭曲。如果观测曲线与模板曲线差异很大，可以适当降低lam；如果 warping 形状过于奇怪，可以提高lam。

4.5 输出与检查

训练结束后，可以检查w的端点和单调性。

print("w[0] =", w[0].item()) print("w[-1] =", w[-1].item()) print("min diff =", torch.min(torch.diff(w)).item())

预期看到w[0]接近 0，w[-1]接近 1，min diff大于 0。由于数值原因，w[0]可能不是严格 0，但应该非常小。

在理想情况下，学习到的w(t)会在t≈0.6时接近 0.35，在t≈0.95时接近 0.70，相当于把观测曲线的时间轴拉回模板时间轴。由于端点约束，w(t)在 0 到 1 之间形成一条单调递增的 S 形曲线。

5. 运行验证：重构误差、对齐误差和可视化

5.1 张量分解的验证方法

张量分解最直接的验证指标是相对重构误差：

mse = torch.mean((cp_reconstruct(A, B, C) - X) ** 2) x_norm = torch.mean(X ** 2) relative_error = torch.sqrt(mse / x_norm) print(f"relative error = {relative_error.item():.4f}")

由于噪声很小，相对误差应该在一个较低水平，例如 0.05 到 0.15。具体数值取决于噪声大小、张量规模和分解秩R。

还可以验证因子矩阵的列和是否为 1：

print(A.sum(dim=0)) print(B.sum(dim=0)) print(C.sum(dim=0))

每一列应该非常接近 1，误差来自 softmax 的数值精度和梯度优化后的余量。如果列和不是 1，说明代码中可能在某个环节破坏了概率约束。

5.2 函数配准的验证方法

函数配准需要同时查看两个指标：对齐误差和 warping 正则度。

align_error = torch.mean((template(w) - g) ** 2).item() smooth_error = torch.mean((w_prime - 1.0) ** 2).item() print(f"align error = {align_error:.6f}") print(f"smooth error = {smooth_error:.6f}")

对齐误差下降说明曲线形态被匹配；平滑误差不能太高，否则说明 warping 过度弯曲。更好的验证方式是把三条曲线画在同一张图上：原始模板f、观测g、配准后的f(w)。配准后的曲线应当和观测曲线基本重叠。

5.3 一个表格记录关键数值

为了快速判断实验是否成功，可以把两个场景的关键指标整理成下面的形式。

运行这些最小案例时，建议先固定随机种子，以便复现。调参时每次只改一个变量，先看损失曲线是否下降，再做下一步。

6. 常见问题排查：从现象到根因

6.1 张量分解相关故障

6.2 函数配准相关故障

6.3 排查链路顺序

遇到问题时，建议按下面的顺序排查，不要直接调大模型或换网络结构。

1. 检查输入数据是否合法，是否包含 NaN 或无穷大。
2. 检查参数化变换是否保持约束，打印中间结果的取值范围。
3. 检查损失函数是否计算正确，分别打印各项损失。
4. 检查优化器参数，尤其是学习率和梯度是否出现异常。
5. 检查初始化和随机种子，确认问题是否在不同初始化下都出现。
6. 最后再检查模型容量、基函数数量和正则权重。

这个顺序能避免很多“找不到原因”的调试过程。大多数概率张量分解和函数配准的问题，根源都在参数化是否光滑、损失权重是否失衡、初始化是否合理这三件事上。

7. 学习环境与生产环境差异

7.1 学习环境：跑通两条最小链路

学习阶段只需要跑通两个最小案例。张量规模可以控制在I=10, J=12, K=15，函数配准的采样点数可以控制在 128。这样的规模在 CPU 上几十秒内就能完成训练，方便观察每一步的中间结果。

学习环境允许更大的容错：可以打开torch.autograd.set_detect_anomaly(True)来定位梯度爆炸的位置，可以在损失中临时去掉正则项观察重构能力，也可以打印每个中间矩阵的形状和数值范围。

7.2 生产环境：还需要额外考虑的问题

一旦进入生产环境，问题会复杂很多。

第一是高维张量的规模。真实张量可能达到1000 × 2000 × 3000，直接存储整个张量并计算einsum是不现实的。这时需要考虑分解算法的在线版本、采样版本或稀疏张量格式，而不是简单扩大维度。

第二是收敛诊断和自动调参。生产环境不能只靠人工看损失曲线，需要设定指标阈值、记录训练日志、监控重构误差和因子退化程度。如果因子矩阵的某个列过于集中，几乎是单点分布，说明模型退化，需要触发告警。

第三是随机种子和可复现性。贝叶斯推断和随机梯度优化都依赖随机性，生产环境必须固定种子或记录种子，否则排查线上问题时无法复现模型。

第四是模型保存和部署。CLR 坐标是优化变量，但部署时通常需要输出概率因子矩阵或 warping 函数本身。因此保存模型时，要同时保存Z和映射函数，或者保存最终概率因子矩阵，并明确标注格式。

下表总结了学习环境与生产环境的重点差异。

8. 最佳实践与扩展方向

8.1 可复用检查清单

下面这份清单可以直接用于任何需要“在单纯形乘积空间上做优化”的项目。

• [ ] 明确约束形式：是单纯形、开单纯形还是带边界的单纯形。
• [ ] 选择合适的映射：softmax 适合生成，CLR 适合优化唯一坐标。
• [ ] 初始化时避免全零，使用小尺度随机初始化。
• [ ] 对输入数据和目标值做标准化或范围检查。
• [ ] 训练过程中定期打印约束是否被满足，比如列和是否为 1。
• [ ] 损失函数至少包含重构项和正则项，且权重可配置。
• [ ] 使用梯度裁剪或 clamp 防止边界附近数值爆炸。
• [ ] 固定随机种子，记录训练参数和版本。
• [ ] 部署前检查最终输出是否仍然满足约束。
• [ ] 生产环境增加日志、监控和异常告警。

8.2 扩展方向

这个技术方向可以延伸到多个更有挑战性的场景。

第一是 Riemannian 优化。CLR 重参数化是一种“坐标变换”，而 Riemannian 优化直接在流形上定义对数映射和指数映射，选择更丰富。如果约束空间不是单纯形而是更一般的流形，比如 Stiefel 流形、Grassmann 流形，那么 Riemannian 工具会更合适。

第二是变分推断中的非欧空间建模。概率张量分解如果引入真正的贝叶斯先验，而不是只在优化中加入一个小 KL 正则项，那么还需要实现完整的变分推断。将 Logistic-Normal 与 CLR 结合，可以导出低方差的梯度估计器。

第三是更复杂的函数配准。实际配准任务中，曲线不一定只有一个波峰，可能还有时间缩放、振幅变化和噪声。可以把 warping 函数与振幅模型联合估计，或者在损失中加入二阶平滑正则，限制 warping 的曲率。

第四是张量分解与函数配准的结合。例如将多条函数型数据组成一个三阶张量，其中包含“样本 × 时间点 × 通道”三个维度。先对时间维度做配准，再对样本和通道维度做概率张量分解，可以用同一个几何框架统一推断。

8.3 练习建议

如果想把这一套方法真正变成自己的工具，建议按顺序完成三个练习。

第一个练习：把张量分解案例中的R从 4 改成 2，观察重构误差如何变化，再改成 8，观察是否过拟合。这个练习能帮助理解分解秩的作用。

第二个练习：在函数配准案例中，把观测曲线的偏移从“整体平移”改成“局部拉伸”，比如第一个峰位置不变，第二个峰向右偏移。观察 warping 函数形状如何变化。

第三个练习：把上面的 CLR 重参数化改写成一个独立的 Python 模块，输入任意矩阵，输出满足概率约束且可微的矩阵结果。这个模块可以直接复用到未来其他带约束的模型里。

这三个练习做完，基本就能掌握“约束空间问题转无约束空间求解”的核心思路。之后再看概率张量分解和函数配准的相关论文，关注点会从“公式是什么”转向“这里的重参数化选择是否合理、数值上是否稳定”，这正是这块内容最有价值的地方。