GLinSAT: The General Linear Satisfiability Neural Network Layer By Accelerated Gradient Descent (NeurIPS, 2024)

Citation

Hongtai Zeng, Chao Yang, Yanzhen Zhou, Cheng Yang, and Qinglai Guo, “GLinSAT: The General Linear Satisfiability Neural Network Layer By Accelerated Gradient Descent,” Advances in Neural Information Processing Systems (NeurIPS), 2024.

Paper · OpenReview · Code

Abstract

GLinSAT is a differentiable neural-network layer that maps model outputs to solutions satisfying bounded general linear constraints. It reformulates the projection as an entropy-regularized linear program, converts the constrained primal problem into a smooth unconstrained dual problem, and solves it using an adaptive accelerated gradient method without matrix factorization. Gradients can be obtained either by differentiating through the iterations or by implicit differentiation with a conjugate-gradient solve. The method is evaluated on constrained traveling-salesman, graph-matching, portfolio-allocation, and unit-commitment tasks.

Introduction

GLinSAT 是一个可以将任意神经网络输出 $c \in R^{n^{'}}$ 映射到满足线性约束的可行解 $x \in R^{n^{'}}$ 的可微映射. 其核心是将原始的线性规划问题 (LP) 转化为一个对偶问题, 并使用 Nesterov 加速梯度下降法求解对偶问题, 从而得到原始问题的最优解.

GLinSAT maps neural network outputs to solutions satisfying general linear constraints

Figure: GLinSAT as a differentiable feasibility layer between a neural network and the downstream loss.

Methodology

Reformulation of the neural network output projection problem

首先考虑一般的约束问题.

对于神经网络输出 $c^{'} \in R^{n^{'}}$ , 我们希望通过可微的方式将其进行规范化, 使得规范后的结果 $x^{'} = GLinSAT (c^{'}) \in R^{n^{'}}$ 满足如下可能的约束条件:
$(1a) (1b) (1c) (1d) A_{1}^{'} x^{'} A_{2}^{'} x^{'} A_{3}^{'} x^{'} x^{'} \leq b_{1}^{'} \geq b_{2}^{'} = b_{3}^{'} \in [l^{'}, u^{'}]$
其中 $A_{1}^{'} \in R^{m_{1}^{'} \times n^{'}}$ , $A_{2}^{'} \in R^{m_{2}^{'} \times n^{'}}$ , $A_{3}^{'} \in R^{m_{3}^{'} \times n^{'}}$ 分别是线性约束的系数矩阵, $b_{1}^{'} \in R^{m_{1}^{'}}$ , $b_{2}^{'} \in R^{m_{2}^{'}}$ , $b_{3}^{'} \in R^{m_{3}^{'}}$ 分别是线性约束的右端向量, $l^{'}$ 和 $u^{'}$ 分别是变量的下界和上界. 并且假设这个问题是可行的.

通过标准化的处理手段引入 slack variables $s_{1} \in R^{m_{1}^{'}}$ , $s_{2} \in R^{m_{2}^{'}}$ , 可以整合为如下标准形式的线性约束 (故后文默认采用标准形式):

(2a) (2b) Ax x = b \in [0, u]

这里, 作者使用内积 $⟨ c, x ⟩$ 衡量 GLinSAT 的输出 $x$ 与输入 $c$ 的相似性.
注意, 文中同时将这个内积叫作 projection, 但是其有别于欧氏投影 ( $∥ c - x ∥^{2} = ∥ c ∥^{2} + ∥ x ∥^{2} - 2 ⟨ c, x ⟩$ ). 因此内积的核心是让输出 $x$ 尽量遵循 $c$ 给出的各个分量的分数排序偏好, 但不要求输出 $x$ 与输入 $c$ 在数值上尽量接近.

因此理论上, 我们可以直接考虑如下的 LP 问题:

maximize subject to ⟨ c, x ⟩ Ax = b x \in [0, u] (LP)

然而一个比较显式的问题是, 对于 LP 问题, 其最优解 $x^{⋆} (c)$ 并不一定是唯一的, 并且通常取在线性约束的多面体构成的顶点 $V = {v_{1}, v_{2}, \dots, v_{k}}$ 上. 因此, LP 事实上是在执行一个 argmax 操作 $x^{⋆} (c) \in ar g v_{k} \in V max ⟨ c, v_{k} ⟩$ 即一个分段线性函数. 故这样的映射 $c \mapsto x^{⋆} (c)$ 是一个几乎处处梯度为 0, 在分段点处不可微的函数, 这对于神经网络的反向传播是一个很大的障碍.

因此, 作者在这里提出了 entropy regularization 进行光滑化处理, 得到的正则后的目标函数为:

minimize_{0 \leq x \leq u} subject to - ⟨ c, x ⟩ + \frac{1}{θ} j = 1 \sum n [\frac{x _{j}}{u _{j}} lo g \frac{x _{j}}{u _{j}} + (1 - \frac{x _{j}}{u _{j}}) lo g (1 - \frac{x _{j}}{u _{j}})] Ax = b (LP-Ent)

对于这个目标函数的理解如下.
- 由于 $0 \leq x_{j} \leq u_{j}$ , 因此得到归一化变量 $p_{j} = x_{j} / u_{j} \in [0, 1]$ .
- 对应的标量函数 (延拓定义 $ϕ (0) = 0 lo g 0 := 0$ ) (注意到其恰为 Bernoulli 的负 entropy)
  $ϕ (p) = p lo g p + (1 - p) lo g (1 - p), p \in [0, 1]$
  其二阶导 $ϕ^{''} (p) = \frac{1}{p ( 1 - p )} \geq 4$ , 因此 $ϕ (p)$ 是一个严格凸函数, 且在 $p = 1/2$ 处取得最小值 $ϕ (1/2) = - lo g 2$ . 其图象如下所示.
  
  Figure: The Bernoulli negative entropy regularizer favors smooth interior solutions. 则正则项可以表达为:
  $\frac{1}{θ} j = 1 \sum n [\frac{x _{j}}{u _{j}} lo g \frac{x _{j}}{u _{j}} + (1 - \frac{x _{j}}{u _{j}}) lo g (1 - \frac{x _{j}}{u _{j}})] = \frac{1}{θ} j = 1 \sum n ϕ (p_{j})$
  也由于 $ϕ$ 的特性, 这个正则项会鼓励每个 $p_{j}$ 尽量接近 $1/2$ , 也就是鼓励每个 $x_{j}$ 尽量接近 $u_{j} /2$ , 即将 $x$ 推向内部的平滑解.
- 若只考虑一维的情况, 则目标函数为
  $0 \leq x \leq u min - c x + \frac{1}{θ} [\frac{x}{u} lo g \frac{x}{u} + (1 - \frac{x}{u}) lo g (1 - \frac{x}{u})] .$
  可以给出对应的最优 closed form solution:
  $x_{θ}^{⋆} (c) = \frac{u}{1 + exp ( - θ u c )} := u σ (θ u c)$
  其中 $σ (\cdot)$ 是 sigmoid 函数. 另外, 观察到 $θ$ 的作用类似于 temperature (严格说是温度的倒数), 控制了得到的输出的平滑以及均匀程度,具体而言:
  - 当 $θ \to \infty$ , 则 $x_{θ}^{⋆} (c) \to u 1_{{c > 0}}$ , 即恢复到原始的 argmax 的非光滑解. 故 $θ$ 越大, 则越接近原始的 LP 解, 但是梯度也越不稳定. (注意, 当 LP 非退化时, 该逼近方才成立)
  - 当 $θ \to 0$ , 则 $x_{θ}^{⋆} (c) \to u /2$ , 由 entropy 直接得到区间中间的正则解. 故 $θ$ 越小, 则越接近平滑解 (越平均), 但是偏离原始的 LP 解也越远.

接着, 我们先暂时只考虑 $LP-Ent$ 的等式约束, 将不等式 box-constrain 暂时作为定义域的限制. 逐步考虑其对偶问题.

引入 Lagrangian Multiplier $y \in R^{m}$ , 则对应的 Lagrangian 为:
$L (x, y) = - ⟨ c, x ⟩ + \frac{1}{θ} j = 1 \sum n ϕ (\frac{x _{j}}{u _{j}}) - ⟨ y, Ax - b ⟩ = - ⟨ c, x ⟩ + \frac{1}{θ} 1^{⊤} [x ⊘ u ⊙ lo g (x ⊘ u) + (1 - x ⊘ u) ⊙ lo g (1 - x ⊘ u)] - ⟨ y, Ax - b ⟩$
- 其中 $⊘$ 表示逐元素除法, $⊙$ 表示逐元素乘法, $1$ 表示全 1 向量.
对应的对偶函数为:
$g (y) = 0 \leq x \leq u in f L (x, y) = 0 \leq x \leq u in f (- ⟨ c, x ⟩ + \frac{1}{θ} j = 1 \sum n ϕ (\frac{x _{j}}{u _{j}}) - ⟨ y, Ax ⟩) + ⟨ y, b ⟩ = 0 \leq x \leq u in f (- ⟨ A^{⊤} y + c, x ⟩ + \frac{1}{θ} j = 1 \sum n ϕ (\frac{x _{j}}{u _{j}})) + ⟨ y, b ⟩$
- 其内层最小化问题
  $0 \leq x \leq u in f (- j = 1 \sum n (A^{⊤} y + c)_{j} x_{j} + \frac{1}{θ} j = 1 \sum n ϕ (\frac{x _{j}}{u _{j}})) := 0 \leq x \leq u in f ψ (x)$
  - 注意到其全部都是逐分量求和的, 且分量间彼此解耦, 故可以独立对每个分量求最小值. 故考虑标量函数
    $ψ_{j} (x_{j}; y) = - (A^{⊤} y + c)_{j} x_{j} + \frac{1}{θ} ϕ (\frac{x _{j}}{u _{j}})$
    其一阶导数为
    $\frac{\partial ψ _{j} ( x _{j} ; y )}{\partial x _{j}} = - (A^{⊤} y + c)_{j} + \frac{1}{θ u _{j}} lo g \frac{x _{j}}{u _{j} - x _{j}} .$
    且有, $lim_{x_{j} \to 0^{+}} ψ_{j}^{'} (x_{j}; y) \to - \infty$ , $lim_{x_{j} \to u_{j}^{-}} ψ_{j}^{'} (x_{j}; y) \to \infty$ . 又根据其二阶导数
    $\frac{\partial ^{2} ψ _{j} ( x _{j} ; y )}{\partial x _{j}^{2}} = \frac{1}{θ} \frac{1}{x _{j} ( u _{j} - x _{j} )} > 0, \forall0 < x_{j} < u_{j} .$
    由此可知, $ψ_{j} (x_{j}; y)$ 是严格凸函数, 且在区间 $(0, u_{j})$ 内部有唯一的最小值, 且对全部 $j$ 都成立.
  - 对于每个分量的最小值, 由一阶导数为零的条件可得
    $- (A^{⊤} y + c)_{j} + \frac{1}{θ u _{j}} lo g \frac{x _{j}}{u _{j} - x _{j}} ⟹ x_{j} (y) = \frac{u _{j}}{1 + exp ( - θ u _{j} ( A ^{⊤} y + c ) _{j} )} = u_{j} \cdot σ (θ u_{j} (A^{⊤} y + c)_{j}) .$
    故用向量整合, 则为:
    $x (y) = u ⊙ σ (θ u ⊙ (A^{⊤} y + c)) .$
- 在成功求出内层最小化问题 $x (y)$ 后, 则可以将其代入对偶函数 $g (y)$ 中, 得到
  $g (y) = \frac{1}{θ} 1^{⊤} lo g σ (θ u ⊙ (- A^{⊤} y - c)) + ⟨ y, b ⟩$
  因此, 对偶问题可以写为:
  $y \in R^{m} min - g (y) = y \in R^{m} min - \frac{1}{θ} 1^{⊤} lo g σ (θ u ⊙ (- A^{⊤} y - c)) - ⟨ y, b ⟩$
  又根据原函数和对偶函数的关系, 可以证明 $- g (y)$ 是 Lipschitz 光滑的, 因此所有梯度下降法都可以用来求解这个对偶问题.
  - 这里强调, 对于内层问题求出的解 $x (y)$ , 其是 $x (y) \in ar g 0 \leq x \leq u min L (x, y)$ 的最小值. 其求解过程保证了 box-constrain $0 \leq x \leq u$ 是一定可以被满足的. 然而对于等式约束 $Ax = b$ , 其仍然是在 Lagrangian 中被当做惩罚项处理的, 因此 $x (y)$ 并不一定满足等式约束 $Ax = b$ .
    - 只有当 $y$ 取到最优解 $y^{⋆}$ 时, 才能保证 $x (y^{⋆})$ 满足等式约束 $Ax = b$ . 理由如下:
      - 原问题的 Lagrangian 为: $L (x, y) = - ⟨ c, x ⟩ + \frac{1}{θ} \sum_{j = 1}^{n} ϕ (\frac{x _{j}}{u _{j}}) - ⟨ y, Ax - b ⟩$ . 对偶问题为 $g (y) = in f_{0 \leq x \leq u} L (x, y) = L (x (y), y)$ . 由对偶问题的最优性条件可知, $\nabla g (y) = \nabla_{y} L (x (y), y) + [\nabla_{x} L (x (y), y)]^{⊤} \frac{\partial x ( y )}{\partial y}$ 而又根据 $\nabla_{x} L (x (y), y) = 0$ 的最优性条件, 则有 $\nabla g (y) = \nabla_{y} L (x (y), y) = - Ax (y) + b .$ 因此对于最优解 $y^{⋆}$ , 有 $\nabla g (y^{⋆}) = 0 ⟹ Ax (y^{⋆}) = b$ , 即 $x (y^{⋆})$ 满足等式约束 $Ax = b$ .

Forward pass in GLinSAT

由上面的推导, 我们得到对于神经网络的输出 $c \in R^{n}$ , GLinSAT 需要向前传播得到满足约束的 $x \in R^{n}$ . 而该约束已经转化成了一个对偶问题, 其形式为:

y \in R^{m} min F (y)

其中 $F (y) = - g (y)$ 是一个 Lipschitz 光滑函数, 且其梯度为

\nabla F (y) = Ax (y) - b, x (y) = u ⊙ σ (θ u ⊙ (A^{⊤} y + c))

故笼统地讲, GLinSAT 的前向传播过程即为:

c \mapsto y^{⋆} \mapsto x (y^{⋆})

故下面的核心是如何高效地求解对偶问题 $min_{y \in R^{m}} F (y)$ .

已知 $F (y)$ 是一个 Lipschitz 光滑函数, 且其梯度为 $\nabla F (y) = Ax (y) - b$ , 其中 $x (y) = u ⊙ σ (θ u ⊙ (A^{⊤} y + c))$ . 故对任意 $y_{1}, y_{2} \in R^{m}$ , 有
$∥\nabla F (y_{1}) - \nabla F (y_{2}) ∥ = ∥ Ax (y_{1}) - Ax (y_{2}) ∥ \leq ∥ A ∥ \cdot ∥ x (y_{1}) - x (y_{2}) ∥ = ∥ A ∥ \cdot ∥ u ⊙ σ (θ u ⊙ (A^{⊤} y_{1} + c)) - u ⊙ σ (θ u ⊙ (A^{⊤} y_{2} + c)) ∥ \leq \frac{θ}{4} j max u_{j}^{2} \cdot ∥ A ∥^{2} \cdot ∥ y_{1} - y_{2} ∥$
- 其中, 最后一个不等式是根据 sigmoid 的性质: $∣ σ (a) - σ (b) ∣ \leq \frac{1}{4} ∣ a - b ∣$ 得到的. 故 $F (y)$ 是一个 Lipschitz 光滑函数, 且其 Lipschitz 常数为 $L = \frac{θ}{4} max_{j} u_{j}^{2} \cdot ∥ A ∥^{2}$ .
给定 $F$ 的 convex & Lipschitz 光滑性质, 因此可以使用 Nesterov 型的 primal-dual 梯度加速方法 (AdaptivePrimalDualAcceleratedGradientDescent) 对这个目标函数进行优化, 以得到 $O (1/ k^{2})$ 的收敛率.
- 之所以不使用 vanilla NAG, 是因为我们的本质期待还是高效的求解 primal 问题 $x^{⋆}$ . 因此我们在对 $F (y)$ 进行优化的过程中, 同时要维护一个 $x^{k}$ 的优化序列, 使得对于对偶问题的加速求解能够有效转化为 primal 问题的加速求解. 此外 Adaptive 这里是说我们在每次迭代中, 通过某种 linesearch 方法, 得到一个更为适配的符合当前迭代函数的 Lipschitz 常数 (记为 $M^{k}$ ), 以便更快地收敛 (而不是使用全局保守的 Lipschitz 下界).
- 下面具体讨论一下 GLinSAT 的 Nesterov 加速算法. 其主体同样采用三阶段的 Nesterov 加速方法:
  $λ^{k + 1} ζ^{k + 1} η^{k + 1} = (1 - τ^{k + 1}) η^{k} + τ^{k + 1} ζ^{k} = ζ^{k} - α^{k + 1} \nabla F (λ^{k + 1}) = (1 - τ^{k + 1}) η^{k} + τ^{k + 1} ζ^{k + 1}$
  - 其中, $η^{k} \in R^{m}$ 是对偶问题的加速序列, $ζ^{k} \in R^{m}$ 是对偶问题的梯度下降序列, $λ^{k + 1} \in R^{m}$ 是进行梯度计算的查询点. 其都是在对偶空间中完成的. 真正的 primal 变量是通过 $x (λ^{k + 1})$ 来计算的.
  - $τ^{k + 1}$ 是 Nesterov 的加速参数, $α^{k + 1}$ 是步长参数, 其可以通过某种 linesearch 方法自适应地选择.
  - 初始化时, 令 $k = 0$ , $η^{0} = ζ^{0} = y^{0}$ . 且 $x^{0} = x (y^{0})$ . 以及 Lipschitz 常数 $M^{0} = L^{0}$ . 额外两个加速过程的迭代参数, $β^{0} = α^{0} = 0$ .
  - 算法最终的停止条件为
    $∥ A x^{k} - b ∥ \leq ϵ$
  - 首先, 计算 $α^{k + 1}$ , 其为如下二次方程 $M^{k} (α^{k + 1})^{2} - α^{k + 1} - β^{k} = 0$ 的正根, 即
    $α^{k + 1} := \frac{1 + 1 + 4 M ^{k} β ^{k}}{2 M ^{k}} .$
    并接着定义
    $β^{k + 1} := β^{k} + α^{k + 1} .$
    通过上述的定义方法, 恒有如下等式成立, 其是 Nesterov 加速方法推导的核心:
    $M^{k} (α^{k + 1})^{2} \equiv β^{k + 1}$
    该表达式在迭代中始终被满足.
  - 接着, 计算加权系数 $τ^{k + 1}$ , 其为
    $τ^{k + 1} := \frac{α ^{k + 1}}{β ^{k + 1}} .$
  - 然后, 计算查询点 $λ^{k + 1}$ , 其为
    $λ^{k + 1} := (1 - τ^{k + 1}) η^{k} + τ^{k + 1} ζ^{k} .$
  - 然后计算查询点对应的 primal 变量 $x (λ^{k + 1})$ , 其为
    $x (λ^{k + 1}) := u ⊙ σ (- θ u ⊙ (- A^{⊤} λ^{k + 1} - c)) .$
    以及梯度信息
    $\nabla F (λ^{k + 1}) = Ax (λ^{k + 1}) - b .$
  - 然后更新对偶变量 $ζ^{k + 1}$ , 其为
    $ζ^{k + 1} := ζ^{k} - α^{k + 1} \nabla F (λ^{k + 1}) .$
  - 然后更新对偶变量 $η^{k + 1}$ , 其为
    $η^{k + 1} := (1 - τ^{k + 1}) η^{k} + τ^{k + 1} ζ^{k + 1} .$
    若将上述更新过程进行合并整理, 可以得到
    $η^{k + 1} = λ^{k + 1} - \frac{1}{M ^{k}} \nabla F (λ^{k + 1}) .$
  - 最后, 进行自适应的尝试更新 Lipschitz 常数 $M^{k + 1}$ . 观察上面的等式, $1/ M^{k}$ 是更新的步长. 如果 $M^{k}$ 太大, 则步长太小, 更新过于保守. 因此我们需要随着迭代的不断进行, 不断尝试得到一个更小更为准确的曲率系数 $M^{k}$ , 以得到更为合理的步长. 事实上, 我们进行尝试的选择依据就是 $F (y) \leq F (x) + ⟨ \nabla F (x), y - x ⟩ + \frac{M}{2} ∥ y - x ∥_{2}^{2}$ . 将具体的迭代点代入, 最终的 linesearch 的条件可以化简整理为:
    $F (η^{k + 1}) \leq F (λ^{k + 1}) - \frac{1}{2 M ^{k}} ∥\nabla F (λ^{k + 1}) ∥^{2}$
    也就是我们本身是希望 $M^{k}$ 选择的越小越好, 但是如果当前步骤不满足上述条件, 则说明 $M^{k}$ 选择的过小, 则将该 $M^{k}$ double, 直到满足上述条件为止. 然后再下一轮的时候再尝试选择一个更小的 $M^{k + 1}$ (取 half), 以便更快地收敛.
    - 当然在实践中还会有一些工程参数, 例如额外引入 tolerance 阈值 $δ > 0$ . 当 $F (η^{k + 1}) > F (λ^{k + 1}) - \frac{1}{2 M ^{k}} ∥\nabla F (λ^{k + 1}) ∥^{2} + δ$ 时, 说明当前 $M^{k}$ 仍然过小, 应拒绝该次尝试并增大 $M^{k}$ . 当充分下降条件连续成功两次时, 才正式尝试更新 $M^{k + 1} = M^{k} /2$ . 这样可以避免过于频繁地更新 $M^{k}$ 而导致不稳定.

Adaptive primal-dual accelerated gradient descent algorithm used in the GLinSAT forward pass

Figure: APDAGD updates the dual variables while maintaining an averaged primal solution.

上述的迭代复杂度为

O (∥ A ∥_{2} max (u) \frac{θR}{ε}),

$θ$ 越大, 意味着 entropy regularization 越弱, 输出更接近原始 LP 的极点解; 但同时对偶目标的光滑常数变大, forward pass 更难求解.

Backward pass in GLinSAT

整体而言 GLinSAT 在模型 pipeline 中的作用如下.

DATA Neural Network c GLinSAT x Loss Function ℓ (x)

若将 GLinSAT 内部再进一步拆开, 则可以得到如下的流程图:

DATA Neural Network c Dual Problem y^{⋆} (c) Primal Solution x (y^{⋆} (c)) Loss Function ℓ (x)

因此, 在 backward pass 中, 我们需要计算 $\frac{\partial ℓ}{\partial c}$ .

AutoGrad

一种简单的方法就是自动调用 PyTorch 等深度学习框架进行自动的反向传播, 由于注意到在我们求解 GLinSAT 内部时, 要进行 $K$ 次的 APDAGD 迭代, 尽管这个迭代过程本身全部是可微的, 但若 unrolling 展开的话会极大的增加计算图的复杂度, 使得反向传播的计算量和显存消耗都非常大.

CG-based Backward

更为经济的做法是探究最终最优点处的模型结构, 利用模型的相关数学方程结构来帮助推导.

回顾, 优化的对偶目标为
$y \in R^{m} min F (y)$
且
$\nabla F (y) = Ax (y) - b$
在最优点处, 根据 $\nabla F (y^{⋆}) = 0$ , 则有
$h (y, c) := Ax (y, c) - b = A (u ⊙ σ (- θ u ⊙ (- A^{⊤} y - c))) - b,$
并且在最优点处满足 $h (y^{⋆}, c) = 0$ .
因此, 在最优点处, 我们有一个隐函数关系 $h (y^{⋆}; c) = 0$ . 根据隐函数定理, 我们可以对 $y^{⋆}$ 关于 $c$ 求导, 得到
$\frac{\partial y}{\partial c} = - (\frac{\partial h}{\partial y})^{- 1} \frac{\partial h}{\partial c} .$
其中可以分别代数求得
$\frac{\partial h}{\partial y} = A Diag (θ x ⊙ (u - x)) A^{⊤} .$
以及
$\frac{\partial h}{\partial c} = A Diag (θ x ⊙ (u - x)) .$
因此最终得到的
$\frac{\partial y}{\partial c} = - [A Diag (θ x ⊙ (u - x)) A^{⊤}]^{- 1} A Diag (θ x ⊙ (u - x)) .$
因此再代入更完整的计算图中, 从 $c$ 出发, 最终得到的损失函数的梯度关系为:
$\frac{d l}{d c} = \frac{\partial l}{\partial x} (\frac{\partial x}{\partial c} + \frac{\partial x}{\partial y} \frac{\partial y}{\partial c}) .$
故代入隐式微分结果, 有
$\frac{\partial l}{\partial c} = \frac{\partial l}{\partial x} \frac{\partial x}{\partial c} - \frac{\partial l}{\partial x} \frac{\partial x}{\partial y} (\frac{\partial h}{\partial y})^{- 1} \frac{\partial h}{\partial c} .$
不过在论文中还进行了一个更 generalized 的处理. 其若允许对偶变量 $y$ 也和 $x$ 一样, 直接参与到损失函数中, 则最终的梯度公式为:
$\frac{\partial l}{\partial c} = \frac{\partial l}{\partial x} \frac{\partial x}{\partial c} - (\frac{\partial l}{\partial x} \frac{\partial x}{\partial y} + \frac{\partial l}{\partial y}) (\frac{\partial h}{\partial y})^{- 1} \frac{\partial h}{\partial c} .$
若想还原回原始的 $y$ 不参与损失函数的情况, 则只需要令 $\frac{\partial l}{\partial y} = 0$ 即可.
在得到理论的梯度传播公式后, 我们还需要考虑如何高效地计算这个 VJP. 故我们可以将其转化为一个线性方程组求解问题, 其形式为
$(\frac{\partial l}{\partial x} \frac{\partial x}{\partial y} + \frac{\partial l}{\partial y}) (\frac{\partial h}{\partial y})^{- 1} ⟹ (\frac{\partial h}{\partial y}) v = (\frac{\partial l}{\partial x} \frac{\partial x}{\partial y} + \frac{\partial l}{\partial y})^{⊤} =: r, Hv = r$
从上述线性系统中求解出 $v$ . 其中 $H := \frac{\partial h}{\partial y}$ , 而列向量 $r$ 已在上式定义.
- 由于 $h (y) = Ax (y) - b$ , 故将具体的表达式代入整理, 则上述线性系统 LHS 的系数为
  $H = \frac{\partial h}{\partial y} = A Diag (θ x ⊙ (u - x)) A^{⊤} := AD A^{⊤} .$
- 且 $\frac{\partial h}{\partial y}$ 是对称半正定的. 若进一步 $A$ 是行满秩, 则是对称正定的. 因此可以使用 CG 方法高效求解该线性系统, 以得到 $v$ . 这是因为, 上面我们展示了, $H = AD A^{⊤}$ 是对称正定的, 因此求解 $Hv = r$ 就等价于求解 QP
  $v min \frac{1}{2} v^{⊤} Hv - r^{⊤} v$
  而这样的方法就很适合使用 CG 方法.
- Conjugate Gradient (CG) 方法是求解对称正定线性系统的经典方法. 其基本思想就是沿着共轭方向去进行搜索迭代, 而共轭方向可以简单理解为由 $H$ 诱导的正交方向. 其具体迭代过程如下:
  - 初始化 $v^{0} = 0$ , 初始化残差 $s^{0} = r - H v^{0} = r$ , 初始化搜索方向 $p^{0} = s^{0}$ . 即第一步沿着负梯度方向进行搜索.
  - 在第 $k$ 步迭代中, 当前迭代点为 $v^{k}$ , 搜索方向为 $p^{k}$ , 步长为 $γ^{k}$ , 我们希望沿着直线 $v^{k + 1} = v^{k} + γ^{k} p^{k}$ 进行搜索, 使得 $v^{k + 1}$ 最小化二次函数 $min_{γ} Q (v) = \frac{1}{2} v^{⊤} Hv - r^{⊤} v$ . 故可以求得第 $k$ 步的步长为 $γ^{k} = \frac{( s ^{k} ) ^{⊤} s ^{k}}{( p ^{k} ) ^{⊤} H p ^{k}} .$ 进而以此更新迭代点 $v^{k + 1} = v^{k} + γ^{k} p^{k} .$ 更新残差 $s^{k + 1} = s^{k} - γ^{k} H p^{k} .$
  - 检查迭代停止条件 $\frac{∥ s ^{k + 1} ∥}{∥ r ∥} \leq ϵ$ 若满足, 则停止迭代, 否则继续迭代, 更新搜索方向: $p^{k + 1} = s^{k + 1} + \frac{( s ^{k + 1} ) ^{⊤} s ^{k + 1}}{( s ^{k} ) ^{⊤} s ^{k}} p^{k} .$
    - 这里 $\frac{( s ^{k + 1} ) ^{⊤} s ^{k + 1}}{( s ^{k} ) ^{⊤} s ^{k}}$ 是 CG 方法中用于构造共轭搜索方向的系数; 在精确算术下, 新方向与此前搜索方向在 $H$ 诱导的内积下保持共轭.
- 可以看到, 通过引入 CG 方法, 我们完全摒弃了矩阵求逆, 而是化成了一系列的矩阵乘法. 而且中间的 $D$ 还是对角矩阵, 其实质上进一步被简化为标量乘法. 计算量被大幅简化.

Experiments

Predictive Portfolio Allocation

对于 predictive portfolio allocation 问题, 模型的目标是最大化 Sharpe ratio, 同时保证输出满足线性的资产配置约束. 其总体 pipeline 如下:

Market Data ⟶ NN ⟶ c ⟶ GLinSAT ⟶ x ⟶ Portfolio Performance

其中 $x \in [0, 1]^{n}$ 是投资组合的权重向量, 并假定如下约束:

资金全部投入, 即 $\sum_{j = 1}^{n} x_{j} = 1$ .
给定偏好资产集合 $S \subseteq {1, 2, \dots, n}$ , 以及配置权重下限 $q \in [0, 1]$ , 要求偏好资产配置份额不低于 $q$ , 即 $\sum_{j \in S} x_{j} \geq q$ . 文中取 $q = 0.5$ .
不允许做空或杠杆, 即 $0 \leq x_{j} \leq 1, \forall j$ .

BYNotes

Explorer