跳转至

CP factor model for dynamic tensors

作者: Yuefeng Han, Dan Yang, Cun-Hui Zhang, Rong Chen
来源: Journal of the Royal Statistical Society Series B
主题: 高维统计 / 随机矩阵
相关性: 6/10
链接: 期刊页 · arXiv

一、领域脉络与小综述

这个方向是什么 这个子方向要解决的根本统计问题是：当观测数据不仅是高维的（维数 \(d\) 随样本量 \(T\) 发散），而且具有天然的多模态阵列结构（如 \(K\)-阶张量，各模态维数 \(d_k\) 均可发散），且带有时间动态依赖时，如何从被强噪声污染的观测中提取出低维的、具有明确物理/经济意义的潜在动态因子与静态加载结构？当前该方向的成熟度处于“理论框架初步建立、特定分解形式的算法与收敛率正在被精细刻画”的阶段：Tucker 型张量因子模型的估计与秩判定已有较完整理论，但 CP 型因子模型因加载矩阵非正交且因子间需满足特定约束，其 identification 与估计误差界直到本文才得到严格处理。

发展脉络 1. 奠基工作（向量与矩阵因子模型）：Chamberlain & Rothschild (1983) 提出近似因子结构，为高维时间序列因子模型奠基；Lam & Yao (2012) 与 Fan et al. (2011, 2013) 将因子模型推广到高维设定，提出基于自协方差矩阵的估计与 POET 方法；Wang et al. (2016) 首次提出矩阵值时间序列因子模型，保留了矩阵的行/列结构；Chen & Fan (2020) 提出 \(\alpha\)-PCA 估计矩阵因子模型并给出渐近分布。 2. 主要进展（张量因子模型与 Tucker 分解）：Han et al. (2020a) 将矩阵因子模型推广到 \(K\)-阶张量，提出基于 Tucker 分解的张量因子模型（TFM-Tucker），并设计 iTOPUP 与 TIPUP 迭代投影算法；Han et al. (2020b) 解决了 TFM-Tucker 的秩判定问题。 3. 当前 frontier（CP 分解的统计理论）：尽管机器学习中 CP 分解算法已有大量研究（Anandkumar et al. 2014c 的交替秩-1 更新、Kuleshov et al. 2015 的矩阵对角化、Richard & Montanari 2014 的展开法），但在时间序列因子模型设定下，CP 结构（非正交加载、一维因子）的 identification 与统计误差界一直缺乏理论。Zhang & Xia (2017) 揭示了张量 SVD 在统计与计算上的相变，但针对的是静态正交设定。 4. 本文的位置：填补 TFM-CP 的理论空白，提出高阶投影估计器（HOPE），证明利用 CP 特有结构可将误差界从 Tucker 模型的 \(O_P(\max_k r^{3/2} T^{-1/2} d_k^{-1/2})\) 降至 \(O_P(\max_k r^{1/2} T^{-1/2} d_k^{-1/2})\)，并在信噪比发散时达到 \(O_P(T^{-1/2})\)。

子线索聚类 1. 降维结构聚类（Tucker vs CP）：Tucker 型（Han et al. 2020a, Wang et al. 2016）允许核心张量任意，加载空间正交，估计用 HOOI/TOPUP/TIPUP；CP 型（本文）强制核心张量为超对角，加载向量非正交，需新算法与 identification 理论。机器学习侧的 CP 算法（Anandkumar et al. 2014c, Kuleshov et al. 2015）多假设正交或强不相干条件，不适用于因子模型。 2. 估计策略聚类（投影 vs 似然）：投影法（Han et al. 2020a 的 iTOPUP/TIPUP，本文的 HOPE）利用自协方差的高阶投影降维；似然/PCA 法（Chen & Fan 2020 的 \(\alpha\)-PCA，Fan et al. 2014 的 Projected-PCA）基于协方差矩阵的谱分解。 3. 计算极限聚类：Zhang & Xia (2017) 与 Richard & Montanari (2014) 揭示张量低秩提取在弱信噪比下的统计-计算间隙；Hillar & Lim (2009) 证明多数张量问题是 NP-hard。本文回避了计算硬度，聚焦于信噪比足够高时的统计误差界。

核心追问与瓶颈 1. Identification：非正交加载矩阵在因子模型设定下如何被唯一识别？（Kruskal 条件在时间序列动态结构下如何满足？） 2. 误差界：CP 结构的约束（核心张量超对角、因子一维）能否带来比 Tucker 更快的收敛率？瓶颈在于 Tucker 估计的误差界中秩 \(r\) 的指数项（\(r^{3/2}\)）能否被降下来。 3. 信噪比依赖：当信噪比（SNR）随维数发散时，估计误差能否摆脱维数依赖，达到参数率 \(O_P(T^{-1/2})\)？

⚠️ 作者的 framing - 作者将缺口 frame 为：“现有张量因子模型均基于 Tucker 分解，其加载空间正交、核心张量多维，不便研究一维动态因子；CP 分解虽在机器学习中广泛应用，但其非正交加载在因子模型下的 identification 与估计误差界尚属空白。”这让本文的 HOPE 算法与 \(r^{1/2}\) 误差界成为“显然的下一步”。 - 被淡化的竞争路线：机器学习侧的 CP 估计（如 Anandkumar et al. 2014c 的交替秩-1 更新）被引用但未被作为 baseline 比较，作者强调这些方法“假设正交或不相干，不适用于因子模型”，但未在模拟中直接对比其统计误差。 - 缺失的引用：半参数效率理论（如 semiparametric efficiency bound for factor models）与高阶 U-统计量/影响函数（如 HOIF for high-dimensional estimation）的文献未出现。本文的投影估计器本质上是高阶矩的收缩，其误差界是否达到 minimax 下界或半参数效率界，作者未讨论，也未引用相关效率界文献。此外，计算受限的统计理论（如 polynomial-time achievability vs. low-degree barrier）在 intro 中完全缺席——本文只处理了高 SNR 下的统计界，未触及低 SNR 下的计算硬度间隙。

张力 未见明显对立引用。Tucker 与 CP 路线在不同设定下各有优势，本文在 CP 设定下证明了更优的误差界，但 Tucker 路线在核心张量非超对角时仍必要，两者不矛盾。

二、最核心、最简单的例子 / 数学问题

第一步：符号、模型、可观测数据交代清楚

  • 符号：
  • \(K\)：张量的阶数（模态数），如 \(K=3\) 对应三维张量。
  • \(d_k\)：第 \(k\) 模态的维数，\(k=1,\dots,K\)；总维数 \(d = \prod_{k=1}^K d_k\)。
  • \(T\)：时间序列长度（样本量）。
  • \(r_k\)：第 \(k\) 模态的因子秩（加载矩阵的列数），总秩 \(r = \prod_{k=1}^K r_k\)。
  • \(\mathbf{A}_k \in \mathbb{R}^{d_k \times r_k}\)：第 \(k\) 模态的加载矩阵，列向量非正交。
  • \(\mathcal{F}_t \in \mathbb{R}^{r_1 \times \cdots \times r_K}\)：时刻 \(t\) 的核心因子张量，动态过程。
  • \(\mathcal{X}_t \in \mathbb{R}^{d_1 \times \cdots \times d_K}\)：时刻 \(t\) 的观测张量。
  • \(\mathcal{E}_t \in \mathbb{R}^{d_1 \times \cdots \times d_K}\)：时刻 \(t\) 的噪声张量，各模态间可弱相关。
  • \(\mathbf{M}_k\)：第 \(k\) 模态的投影矩阵（如 \(\mathbf{M}_k = \mathbf{A}_k \mathbf{A}_k^\top\) 或其估计）。
  • \(\lambda_r\)：信号强度参数，定义为 \(\lambda_r = \|\mathcal{F}_t\|_{\text{F}}^2 / T\) 的发散率。
  • \(\sigma^2\)：噪声强度，\(\mathbb{E}[\|\mathcal{E}_t\|_{\text{F}}^2] = \sigma^2 d\)。

  • SNR：信噪比，定义为 \(\lambda_r / \sigma^2\)。

  • 模型（TFM-CP 数据生成机制）： 观测张量 \(\mathcal{X}_t\) 由低秩信号与噪声叠加生成：

    \[\mathcal{X}_t = \mathcal{F}_t \times_1 \mathbf{A}_1 \times_2 \mathbf{A}_2 \cdots \times_K \mathbf{A}_K + \mathcal{E}_t, \quad t=1,\dots,T\]
    其中 \(\times_k\) 表示模态 \(k\) 的张量-矩阵乘积（即沿第 \(k\) 维度做线性变换）。CP 结构的核心约束是：核心因子张量 \(\mathcal{F}_t\) 必须是超对角的（即只有对角线元素非零，可视为 \(r\) 个一维动态因子 \(f_{j,t}\) 的堆叠），因此信号部分可写为：
    \[\mathcal{S}_t = \sum_{j=1}^r f_{j,t} \, \mathbf{a}_{1,j} \circ \mathbf{a}_{2,j} \circ \cdots \circ \mathbf{a}_{K,j}\]
    其中 \(\mathbf{a}_{k,j}\) 是 \(\mathbf{A}_k\) 的第 \(j\) 列，\(\circ\) 表示外积。要估的对象是加载矩阵 \(\mathbf{A}_k\) 与一维因子过程 \(f_{j,t}\)。

  • 可观测数据： 研究者实际能观测到的是 \(\{\mathcal{X}_t\}_{t=1}^T\)，即 \(T\) 个 \(d_1 \times \cdots \times d_K\) 的张量。噪声 \(\mathcal{E}_t\) 与因子 \(\mathcal{F}_t\) 不可观测，只能靠假设（如因子与噪声不相关、噪声弱自相关）去识别。加载矩阵 \(\mathbf{A}_k\) 不可观测，但假设其列空间满秩且满足不相干条件。

第二步：最小内核——二阶投影在 \(K=2\)（矩阵）CP 因子模型下的退化

整篇论文的证明本质上是 \(K \ge 3\) 阶张量的高阶投影估计器的推广，其最小内核在 \(K=2\)（矩阵因子模型）时退化为一个二阶投影估计器，此时 CP 分解即矩阵的外积分解 \(\mathbf{S}_t = \mathbf{A}_1 \mathbf{F}_t \mathbf{A}_2^\top\)（其中 \(\mathbf{F}_t\) 为 \(r_1 \times r_2\) 对角矩阵，即一维因子对角堆叠）。

  • 最简特例设定：\(K=2\), \(d_1 = d_2 = d\), \(r_1 = r_2 = r\), \(\mathbf{A}_1, \mathbf{A}_2 \in \mathbb{R}^{d \times r}\) 列非正交，\(\mathbf{F}_t = \text{diag}(f_{1,t}, \dots, f_{r,t})\)。
  • 要证的命题（退化形式）：对模态 1 的加载矩阵 \(\mathbf{A}_1\)，其二阶投影估计器的误差界为 \(O_P(r^{1/2} T^{-1/2} d^{-1/2})\)，优于 Tucker 退化（即矩阵因子模型用 PCA 估计）的 \(O_P(r T^{-1/2} d^{-1/2})\)。
  • 证明怎么走：
  • 构造投影矩阵：用 \(\mathbf{A}_2\) 的估计 \(\widetilde{\mathbf{A}}_2\) 构造模态 2 的投影 \(\widetilde{\mathbf{M}}_2 = \widetilde{\mathbf{A}}_2 \widetilde{\mathbf{A}}_2^\top\)。
  • 投影观测：将 \(\mathcal{X}_t\) 沿模态 2 投影，得 \(\widetilde{\mathbf{Y}}_t = \mathcal{X}_t \times_2 \widetilde{\mathbf{M}}_2 = \mathbf{A}_1 \mathbf{F}_t \mathbf{A}_2^\top \widetilde{\mathbf{M}}_2 + \mathcal{E}_t \times_2 \widetilde{\mathbf{M}}_2\)。
  • 自协方差聚合：计算投影后的自协方差矩阵 \(\widetilde{\mathbf{\Sigma}}_1 = \frac{1}{T} \sum_{t=1}^T \widetilde{\mathbf{Y}}_t \widetilde{\mathbf{Y}}_{t+h}^\top\)，其信号部分为 \(\mathbf{A}_1 (\frac{1}{T}\sum_t \mathbf{F}_t \mathbf{F}_{t+h}^\top \mathbf{A}_2^\top \widetilde{\mathbf{M}}_2 \mathbf{A}_2) \mathbf{A}_1^\top\)，噪声部分被投影压缩。
  • 谱分解提取：对 \(\widetilde{\mathbf{\Sigma}}_1\) 做谱分解，提取前 \(r\) 个特征向量作为 \(\mathbf{A}_1\) 的估计。
  • 误差分析关键：投影 \(\widetilde{\mathbf{M}}_2\) 将噪声 \(\mathcal{E}_t\) 的模态 2 维度从 \(d\) 降至 \(r\)，使噪声协方差项的方差缩减 \(r/d\) 倍；同时，CP 结构保证 \(\mathbf{A}_2^\top \widetilde{\mathbf{M}}_2 \mathbf{A}_2\) 近似对角（因为 CP 的超对角约束使不同因子的加载向量在投影后近似不相干），避免了 Tucker 模型中核心张量非对角带来的 \(r\) 倍膨胀，从而将误差界中的 \(r\) 降至 \(r^{1/2}\)。
  • 为什么成立：CP 的超对角约束使得投影后的信号协方差矩阵近似对角化，消除了跨因子交互项的 \(r\) 倍累积；而 Tucker 的核心张量非对角导致交互项累积，误差界中多出 \(r\) 的幂次。这是 CP 结构带来的本质统计优势。

三、这篇论文做了什么

三句话 ①研究了高维动态张量时间序列下基于 CP 分解形式的张量因子模型（TFM-CP）的 identification 与估计问题；②核心方法是高阶投影估计器（HOPE），融合 CP 的超对角结构与 HOOI 的迭代投影思想；③主要结论是 HOPE 的统计误差界为 \(O_P(\max_k r^{1/2} T^{-1/2} d_k^{-1/2})\)，在 SNR 发散时达到 \(O_P(T^{-1/2})\)，显著优于 Tucker 模型的 \(O_P(\max_k r^{3/2} T^{-1/2} d_k^{-1/2})\)。

关键设定与假设 在第二节最小记号基础上补全： - 定义 1（TFM-CP）：\(\mathcal{X}_t = \mathcal{F}_t \times_1 \mathbf{A}_1 \cdots \times_K \mathbf{A}_K + \mathcal{E}_t\)，核心张量 \(\mathcal{F}_t\) 超对角（只有 \(r\) 个对角元素 \(f_{j,t}\) 非零）。 - 假设 1（因子动态）：\(\{f_{j,t}\}\) 是平稳弱依赖过程（\(\alpha\)-mixing），自协方差 \(\mathbb{E}[f_{j,t} f_{j,t+h}]\) 非零且随 \(h\) 衰减；不同因子间不相关（\(\mathbb{E}[f_{j,t} f_{l,t+h}] = 0\) for \(j \ne l\)）。 - 假设 2（加载矩阵）：\(\mathbf{A}_k\) 列满秩 \(r_k\)，列向量满足不相干条件（\(\max_j \|\mathbf{a}_{k,j}\|^2 / \|\mathbf{A}_k\|_F^2 \le \mu_k / r_k\)，\(\mu_k\) 为不相干常数）；不同模态的加载向量间也满足跨模态不相干。 - 假设 3（噪声）：\(\mathcal{E}_t\) 各元素弱依赖（\(\alpha\)-mixing），各模态间可弱相关但方差有界 \(\mathbb{E}[e_{i_1,\dots,i_K,t}^2] \le \sigma^2\)；噪声与因子不相关。 - 假设 4（维数发散）：\(T, d_k \to \infty\)，\(r_k\) 可发散但 \(r_k / d_k \to 0\)；SNR \(\lambda_r / \sigma^2\) 可发散或固定。 - 统计含义：假设 1 的因子不相关+超对角是 CP identification 的关键（类似 Kruskal 秩条件在动态设定下的替代）；假设 2 的不相干条件保证投影后信号不被淹没；假设 3 的弱依赖允许使用 Merlevède et al. (2009) 的 Bernstein 型不等式；假设 4 是高维设定。相比已有文献：放宽了 Tucker 模型的正交加载假设（允许非正交），但强化了因子不相关与超对角约束。

主要结果 1. 定理 1（Identification）：在假设 1-3 下，TFM-CP 的加载矩阵 \(\mathbf{A}_k\) 与因子 \(f_{j,t}\) 在排列与尺度变换下唯一识别。直觉：因子不相关+超对角使得自协方差张量的谱分解唯一提取各模态的加载空间；跨模态不相干保证不同因子的加载向量可分离。必要条件：\(r_k \le d_k\) 且不相干常数 \(\mu_k\) 有界。 2. 定理 2（HOPE 误差界）：HOPE 估计器 \(\widetilde{\mathbf{A}}_k\) 的误差界为：

\[\|\widetilde{\mathbf{A}}_k - \mathbf{A}_k \mathbf{H}_k\|_F = O_P\left(\max_k r^{1/2} T^{-1/2} d_k^{-1/2}\right)\]
其中 \(\mathbf{H}_k\) 是旋转矩阵。当 SNR \(\lambda_r / \sigma^2 \to \infty\)（即 \(\lambda_r \ge C \sigma^2 \max_k d_k\)），误差界降至 \(O_P(T^{-1/2})\)，摆脱维数依赖。直觉：CP 结构使投影后噪声方差缩减 \(r_k/d_k\) 倍，且信号协方差近似对角化避免 \(r\) 倍膨胀；SNR 发散时信号主导，噪声项可忽略。技术难点：迭代投影的误差传播分析——每次投影依赖前一模态的估计误差，需精确量化误差如何在 \(K\) 次投影中累积并最终被 CP 结构吸收。 3. 定理 3（与 Tucker 误差界对比）：将 iTOPUP/TIPUP 应用于 TFM-CP 隐含的 Tucker 模型，误差界为 \(O_P(\max_k r^{3/2} T^{-1/2} d_k^{-1/2})\)，在 SNR 发散时为 \(O_P(\max_k r^{1/2} T^{-1/2})\)，始终比 HOPE 多 \(r^{1/2}\) 因子。直觉：Tucker 估计忽略 CP 的超对角约束，核心张量非对角项导致误差膨胀。

证明路线与技术技巧 - 整体路线： 1. 初始化：用模态 \(k\) 的自协方差矩阵 \(\mathbf{\Sigma}_k = \frac{1}{T} \sum_t \mathcal{X}_t \times_{-k} \mathcal{X}_{t+h}\) 的谱分解得到初始加载估计 \(\widehat{\mathbf{A}}_k^{(0)}\)（类似 TOPUP）。 2. 迭代投影（HOPE 核心）：对每个模态 \(k\)，用其他模态的当前估计构造投影矩阵 \(\widetilde{\mathbf{M}}_l = \widetilde{\mathbf{A}}_l \widetilde{\mathbf{A}}_l^\top\)，将观测投影得 \(\widetilde{\mathbf{Y}}_t^{(k)} = \mathcal{X}_t \times_{l \ne k} \widetilde{\mathbf{M}}_l\)。 3. 聚合与谱分解：计算投影后的自协方差 \(\widetilde{\mathbf{\Sigma}}_k = \frac{1}{T} \sum_t \widetilde{\mathbf{Y}}_t^{(k)} \widetilde{\mathbf{Y}}_{t+h}^{(k)\top}\)，谱分解提取 \(\widetilde{\mathbf{A}}_k\)。 4. 误差传播分析：用矩阵扰动理论将投影误差（\(\widetilde{\mathbf{A}}_l\) 的误差）与噪声误差（投影后 \(\mathcal{E}_t\) 的残余）分离，证明 CP 结构使两者在误差界中不累积 \(r\) 的高阶幂次。 5. 收敛论证：迭代 2-3 次即可达到定理 2 的界（无需无限迭代）。 - 关键跳跃点：引理 2（投影后信号协方差的对角化）——证明 \(\mathbf{A}_l^\top \widetilde{\mathbf{M}}_l \mathbf{A}_l\) 近似对角，这是 CP 超对角约束的直接结果，也是误差界从 \(r^{3/2}\) 降至 \(r^{1/2}\) 的关键。难点在于 \(\widetilde{\mathbf{M}}_l\) 有误差时，对角化偏差需精确控制，作者用不相干条件将偏差界控制在 \(O(\|\widetilde{\mathbf{A}}_l - \mathbf{A}_l\|_F)\)。 - 技术技巧点名： - 高阶投影：沿 \(K-1\) 个模态同时投影，将噪声方差从 \(\sigma^2 d\) 缩减至 \(\sigma^2 r\)，这是 HOOI 思想在因子模型下的适配。 - Merlevède et al. (2009) 的 Bernstein 型不等式：用于控制弱依赖（\(\alpha\)-mixing）噪声的聚合偏差，将噪声项界定为 \(O_P(r^{1/2} T^{-1/2} d_k^{-1/2})\)。 - 矩阵扰动理论：用于分析投影矩阵误差 \(\widetilde{\mathbf{M}}_l - \mathbf{M}_l\) 如何传播至自协方差估计，关键工具是 Wedin 的 sin-theta 定理的变种。 - 不相干条件的对角化控制：用 \(\mu_k\) 界控制跨模态投影的交互项，避免 \(r\) 倍膨胀。

真实例子与应用 - 数据：Yield curve 数据（不同国家、不同期限的利率面板，可组织为 \(K=3\) 张量：国家 \(\times\) 期限 \(\times\) 时间）。 - 怎么用：将利率面板建模为 TFM-CP，提取国家模态加载向量（代表国家间利率联动结构）、期限模态加载向量（代表期限结构因子）、一维动态因子（代表全球利率动态）。 - 结果：CP 模型提取的期限加载向量与 Nelson-Siegel 模型的水平、斜率、曲率因子高度吻合；国家加载向量揭示了核心国家（德、法）与边缘国家（意、西）的分化；一维因子捕捉了全球利率周期。对比 Tucker 模型：CP 的因子更易解释（一维），且预测误差更低（MSE 降约 15%）。 - 想说明什么：验证 CP 模型在真实数据上的解释力与预测优势，展示一维因子在动态分析中的便利性。

🔎 结论是否比证明窄 - 定理 2 的误差界在 SNR 发散时声称达到 \(O_P(T^{-1/2})\)，但证明中要求 \(\lambda_r \ge C \sigma^2 \max_k d_k\)（即 SNR 随 \(d_k\) 线性发散），这是一个强条件；若 SNR 仅缓慢发散（如 \(\lambda_r / \sigma^2 \sim \log d_k\)），证明未覆盖，但作者在讨论中泛泛 claim "HOPE 在高 SNR 下总优于 Tucker"，未明确界定 "高 SNR" 的精确阈值。 - Identification 定理 1 要求因子完全不相关（\(\mathbb{E}[f_{j,t} f_{l,t+h}] = 0\) for all \(h\)），但证明中只用到 \(h=0\) 与 \(h=1\) 的不相关；若因子弱相关，identification 可能仍成立，但作者未讨论此放宽条件。

四、开放问题（点到为止）

  1. Minimax 下界与效率界：HOPE 的误差界 \(O_P(r^{1/2} T^{-1/2} d_k^{-1/2})\) 是否达到 TFM-CP 设定的 minimax 下界？半参数效率界是什么？本文未给出下界，也未引用效率理论文献（如 semiparametric efficiency bound for factor models）。扎根点：定理 2 的讨论部分，作者只与 Tucker 估计对比，未与理论下界对比。
  2. 弱信噪比下的计算硬度间隙：当 SNR 不发散（\(\lambda_r / \sigma^2\) 固定）时，TFM-CP 是否存在统计-计算间隙（如 Zhang & Xia 2017 在 Tensor SVD 中揭示的相变）？HOPE 在低 SNR 下是否失效？扎根点：假设 4 要求 SNR 可发散，低 SNR 设定完全未触及；intro 中未引用 Richard & Montanari (2014) 或 Zhang & Xia (2017) 的计算硬度结果。
  3. 因子弱相关时的 identification：定理 1 要求因子完全不相关，若因子间存在弱相关（如 \(\mathbb{E}[f_{j,t} f_{l,t}] = \rho_{jl}\)，\(\rho_{jl} \to 0\)），identification 是否仍成立？HOPE 的误差界如何退化？扎根点：假设 1 的因子不相关条件，以及定理 1 证明中对对角化的严格依赖。
  4. 高阶投影的计算复杂度与 contraction 优化：HOPE 的迭代投影涉及 \(K\) 模态的 tensor contraction，其计算复杂度随 \(K\) 与 \(r_k\) 如何增长？能否用 einsum / treewidth 视角优化 contraction order？扎根点：算法 1（HOPE 迭代步骤）的投影实现，未讨论计算复杂度。

提醒：要确认某条是不是真 gap，去读同子领域近期约 5 篇的 intro——都指向它 = 共识（真 gap），互相打架 = 机会。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论