libigl第三章：矩阵和线性代数

Libigl 在密集和稀疏线性代数例程中严重依赖 Eigen 库。 除了几何处理例程外，libigl 还具有引导 Eigen 的线性代数例程，使其感觉更类似于 Matlab 等高级代数库。

切片

Matlab 中一个非常熟悉且功能强大的例程是数组切片。 这允许读取或写入可能不连续的子矩阵。 让我们考虑一下 Matlab 代码：

B = A(R,C);

如果A是一个$m\times n$的矩阵，而R是一个长度为j的行索引的向量，C是一个长度为k的列索引的向量，那么B就是一个$j\times k$的矩阵，它的元素是根据R和C中的索引从A中提取出来的。在libigl中，相同的操作通过slice函数给出：

VectorXi R,C;
MatrixXd A,B;
...
igl::slice(A,R,C,B);

请注意，A 和 B 也可以是稀疏矩阵。
同样，考虑 Matlab 代码：

A(R,C) = B;

这段代码的含义是一个$j\times k$的矩阵B将作为一个子矩阵，通过R和C提供的索引写入A中。在libigl中，相同的操作通过slice_into函数给出：

igl::slice_into(B,R,C,A);

排序

Matlab 和其他高级语言使得提取排序和比较例程的索引变得非常容易。 例如在 Matlab 中，可以这样写：

[Y,I] = sort(X,1,'ascend');

所以，如果X是一个$m\times n$的矩阵，那么Y也是一个$m\times n$的矩阵，其中Y是由X每一列元素按照降序排列后的结果。第二个输出记录了X是如何排序生成Y的，它们存在着这样的对应关系$Y(i,j)=X(I(i,j),j)$.I记录的是一个索引矩阵
同样的功能在libigl中也被支持

igl::sort(X,1,true,Y,I);

类似地，可以在 Matlab 中使用以下命令对整行进行排序：

[Y,I] = sortrows(X,'ascend');

其中，I为长度为m的行索引的向量，其中$Y=X(I,:)$,
在 libigl 中，由以下的代码实现

igl::sortrows(X,true,Y,I);

其中，I同样表示的是排序的索引，可以通过igl::slice(X,I,1,Y)，获得相同的Y;
libigl 中提供了类似的函数：max、min 和 unique。

其他的matlab形式的函数

函数名	描述
igl::all	是否矩阵中的所有元素均为非零元素
igl::any	是否矩阵中的任意元素都是非零元素
igl::cat	连接两个矩阵(在处理稀疏矩阵时很有用)
igl::ceil	舍入元素为最接近的整数
igl::cumsum	矩阵元素的累积和
igl::colon	像 Matlab 的 :, 类似于 Eigen 的 LinSpaced
igl::components	连接图的部件(和matlab的graphconncomp比较)
igl::count	对矩阵行或是列中的非零元素进行计数
igl::cross	每一行的叉乘
igl::cumsum	累计和
igl::dot	每一行的点积
igl::eigs	解决稀疏特征值的问题
igl::find	寻找非零元素的下标
igl::floor	对矩阵的元素进行舍入以取整
igl::histc	对构建直方图时的出现次数进行计数
igl::hsv_to_rgb	将hsv表示的颜色转换为rgb表示（类似于matlab的hsv2rgb）
igl::intersect	设置矩阵元素的交集
igl::isdiag	判断矩阵是否为对角矩阵
igl::ismember	判断A中的元素是否出现在B中
igl::jet	沿着彩虹的量化颜色
igl::max	计算每一行或是列的最大的元素
igl::median	计算每一列的中位数
igl::min	计算每一行或是列的最小值
igl::mod	计算每一列的模
igl::null	计算矩阵的零空间基
igl::nchoosek	计算长度为n的向量的所有长度为k的组合
igl::orth	基的正交化
igl::parula	生成从蓝色到黄色的量化颜色图
igl::pinv	计算 Moore-Penrose 伪逆
igl::rgb_to_hsv	将 RGB 颜色转换为 HSV（参见 Matlab 的 rgb2hsv）
igl::repmat	沿列和行重复矩阵
igl::round	每元素四舍五入到整数
igl::setdiff	设置矩阵元素的差异
igl::setunion	设置矩阵元素的并集
igl::setxor	设置矩阵元素的异“或”
igl::slice	使用索引列表对矩阵的部分进行切片：（对比Matlab的 B = A(I,J)）
igl::slice_mask	使用布尔掩码对矩阵的部分进行切片：（对比matlab的B = A(M,N)）
igl::slice_into	使用索引列表对矩阵赋值的左侧进行切片（参见 Matlab 的 B(I,J) = A）
igl::sort	对矩阵的元素或行进行排序
igl::speye	单位初始化稀疏矩阵
igl::sum	沿列或行求和（稀疏矩阵）
igl::unique	提取矩阵的唯一元素或行

拉普拉斯方程(这块内容需要后续深入研究)

几何处理中常见的线性系统是拉普拉斯方程：
$\Delta z=0$
受限于一些边界条件，例如 Dirichlet （狄利克雷）边界条件（固定值）
$z{|}_{\partial S}=z_{bc}$
在离散设置中，线性系统可以写成：
$\mathbf{L}\mathbf{z}=0$
其中$\mathbf{L}$是一个$n\times n$的拉普拉斯算子，$\mathbf{z}$则是由每个顶点值所组成的向量，大多数$\mathbf{z}$对应于内部顶点并且是未知的，但有些$\mathbf{z}$表示边界顶点处的值。它们的值是已知的，因此我们可以将它们的对应项移到右侧。
从概念上讲，如果我们已经将$\mathbf{z}$排序，很容易实现在$\mathbf{z}$中内部顶点首先出现，然后是边界顶点。
$\left(\begin{array}{cc}{\mathbf{L}}_{in,in}& {\mathbf{L}}_{in,b}\\ {\mathbf{L}}_{b,in}& {\mathbf{L}}_{b,b}\end{array}\right)\left(\begin{array}{c}{\mathbf{z}}_{in}\\ {\mathbf{z}}_b\end{array}\right)=\left(\begin{array}{c}{0}_{in}\\ {\mathbf{z}}_{bc}\end{array}\right)$
方程的底部块不再有意义，所以我们只考虑顶部块：
$\left(\begin{array}{cc}{\mathbf{L}}_{in,in}& {\mathbf{L}}_{in,b}\end{array}\right)\left(\begin{array}{c}{\mathbf{z}}_{in}\\ {\mathbf{z}}_b\end{array}\right)=0_{in}$
我们可以将已知值移动到右侧：
${\mathbf{L}}_{in,in}{\mathbf{z}}_{in}={-\mathbf{L}}_{in,b}{\mathbf{z}}_b$
最后，我们可以求解内部顶点处未知值的方程。
但是，我们的顶点通常不会以这种方式排序。一种选择是对 V 进行排序，然后按上述方法进行操作，然后对解 Z 进行排序以匹配 V。但是，该解不是很通用。
使用数组切片不需要显式排序。相反，我们可以切出子矩阵块（例如${\mathbf{L}}_{in,in},{\mathbf{L}}_{in,b}$）并直接遵循上面的线性代数。然后我们可以将解分割成 Z 对应于内部顶点的行

二次能量最小化

相同的拉普拉斯方程可以通过在相同边界条件下最小化狄利克雷能量来等效地导出：
$mininiz{e}_z\frac{1}{2}{\int }_S||\nabla z|{|}^{2}dA$
在我们的离散网格上，回想一下，这变成
$minimiz{e}_z\frac{1}{2}{\mathbf{z}}^T{\mathbf{G}}^T\mathbf{G}\mathbf{z}\to minimiz{e}_z{\mathbf{z}}^T\mathbf{L}\mathbf{z}$
在受固定值边界条件约束的网格上最小化一些能量的一般问题是如此广泛，以至于 libigl 有一个专用的 api 来解决此类系统。
让我们考虑一个受到不同公共约束的一般二次最小化问题：
$minimiz{e}_z\frac{1}{2}{\mathbf{z}}^T\mathbf{Q}\mathbf{z}+{\mathbf{z}}^T\mathbf{B}+constant$
约束于：
${\mathbf{z}}_b={\mathbf{z}}_{bc}and{\mathbf{A}}_{eq}\mathbf{z}={\mathbf{B}}_{eq}$
其中：

• $\mathbf{Q}$通常是一个$n\times n$的稀疏矩阵，为二次系数的半正定矩阵（Hessian）
• $\mathbf{B}$是一个$n\times 1$的线性系数向量
• ${\mathbf{z}}_b$是$\mathbf{z}$的$|b|\times 1$的一部分，对应于边界或是固定的顶点
• ${\mathbf{z}}_{bc}$是对应于${\mathbf{z}}_b$的已知量的$|b|\times 1$的向量
• ${\mathbf{A}}_{eq}$通常是一个$m\times n$的稀疏矩阵，为线性等式约束系数矩阵（每个约束一行）
• ${\mathbf{B}}_{eq}$是线性等式约束右侧值的一个$m\times 1$的向量

因为我们把${\mathbf{A}}_{eq}\mathbf{z}={\mathbf{B}}_{eq}$的行写作${\mathbf{z}}_b={\mathbf{z}}_{bc}$,所以使得这个规定变得过于笼统了，但是这些固定值约束经常出现，以至于它们值得在 API 中专门放置。
在 libigl 中，求解此类二次优化问题分为两个例程：预计算和求解。预计算仅取决于二次系数、已知值索引和线性约束系数：

igl::min_quad_with_fixed_data mqwf;
igl::min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);

输出是一个结构 mqwf，它包含系统矩阵分解，并在求解过程中使用任意线性项、已知值和右侧的约束：

igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);

其中输出V是一个$n\times 1$的解向量，它的固定值被正确的放置以匹配网格的顶点V。

线性等式约束

我们在上面看到 libigl 中的 min_quad_with_fixed_* 提供了一种解决一般二次规划的紧凑方法。让我们考虑另一个例子，这次是主动线性等式约束。
具体来说，让我们求解双拉普拉斯方程或等效地最小化拉普拉斯能量：
${\Delta }^{2}z=0\leftrightarrow minimiz{e}_z\frac{1}{2}{\int }_S(\Delta z{)}^{2}dA$
受固定值约束和线性等式约束：
$z_a=1,z_b=-1and{z}_c=z_d$
请注意，我们可以用上面熟悉的形式重写最后一个约束：
$z_c-z_d=0$
现在我们可以将Aeq组装成一个$1\times n$的稀疏矩阵，它的在d中与顶点c和a-1对应的列的系数1为1。而右手边的Beq就是简单的0。
在内部，min_quad_with_fixed_* 使用拉格朗日乘子法求解。此方法为每个线性约束添加附加变量（通常是变量$\lambda$的一个$m\times 1$的向量），然后，解决鞍点的问题。
$findsaddl{e}_{z,\lambda }\frac{1}{2}{\mathbf{z}}^T\mathbf{Q}\mathbf{z}+{\mathbf{z}}^T\mathbf{B}+constant+{\lambda }^T({\mathbf{A}}_{eq}{\mathbf{z}-\mathbf{B}}_{eq})$
这可以通过将z和$\lambda$堆叠入一个$(m+n)\times 1$的未知数向量中，用更熟悉的形式重写：
$findsaddl{e}_{z,\lambda }\frac{1}{2}（{\mathbf{z}}^T{\lambda }^T)\left(\begin{array}{cc}\mathbf{Q}& {\mathbf{A}}_{eq}^T\\ {\mathbf{A}}_{eq}& 0\end{array}\right)\left(\begin{array}{c}\mathbf{z}\\ \lambda \end{array}\right)+({\mathbf{z}}^T{\lambda }^T)\left(\begin{array}{c}\mathbf{B}\\ {-\mathbf{B}}_{eq}\end{array}\right)+constant$
区别于$（{\mathbf{z}}^T{\lambda }^T)$揭示了一个线性系统，我们可以解决$\mathbf{z}$和$\lambda$。与直接二次最小化系统的唯一区别是，这个鞍问题系统不会是正定的。因此，我们必须使用不同的分解技术（LDLT 而不是 LLT）：libigl 的 min_quad_with_fixed_precompute 在存在线性等式约束时自动选择正确的求解器。

二次规划

通过允许不等式约束，我们可以进一步推广上一节中的二次优化。 特别是包围盒约束（设置下限和上限）：
$\mathbf{l}<\mathbf{z}<\mathbf{u}$
其中$\mathbf{l}$和$\mathbf{u}$是上边界和下边界的$n\times 1$的向量，是一般线性不等式约束：
${\mathbf{A}}_{ieq}\mathbf{z}<{\mathbf{B}}_{ieq}$
其中${\mathbf{A}}_{ieq}$为$k\times n$的线性系数矩阵，${\mathbf{B}}_{ieq}$则是一个$k\times 1$的右侧约束矩阵
同样，我们过于笼统，因为盒子约束可以写成线性不等式约束的行，但边界出现的频率足以值得一个专用的 api。
Libigl 实现了自己的活动集例程来求解二次规划 (QP)。 该算法通过将违反的不等式约束强制执行为等式和“停用”不再需要的约束来迭代地“激活”违反的不等式约束。
在决定了哪些约束在每次迭代中是活动的之后，问题就简化为受线性等式约束的二次最小化，并调用上一节中的方法。 如此重复直到收敛。
目前，该实现对于框约束和稀疏非重叠线性不等式约束是有效的。
与其他内点方法不同，活动集方法受益于热启动（对解向量$\mathbf{z}$的初始猜测）

igl::active_set_params as;
// Z is optional initial guess and output
igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);

特征值分解

Libigl 对提取广义特征值问题的特征对有基本的支持：
$\mathbf{A}x=\lambda \mathbf{B}x$
其中$\mathbf{A}$是一个稀疏对称矩阵，而$\mathbf{B}$则是一个正定矩阵。在几何处理中，一种常见的情形是，我们使得$\mathbf{A}=\mathbf{L}$，即余切拉普拉斯算子，同时令$\mathbf{B}=\mathbf{M}$即每个顶点的质量矩阵。通常应用程序将使用低频特征值模式。
似于傅里叶分解，面上的函数 f 可以通过其对 Laplace-Beltrami 的特征模式的光谱分解来表示：
$f={\sum }_{i=1}^{\infty }{\alpha }_i{\varphi }_i$
其中每个${\varphi }_i$是一个特征函数，满足：$\Delta {\varphi }_i={\lambda }_i{\varphi }_i$,而且${\alpha }_i$是标量系数。
对于离散三角形网格，存在完全类似的分解，尽管总和有限：
$\mathbf{f}={\sum }_{i=1}^n{\alpha }_i{\varphi }_i$
现在$\mathbf{f}\in R^n$是顶点值的列向量,指定分段线性函数，而${\varphi }_i\in R^n$则是满足一下关系的特征值向量：
$\mathbf{L}{\varphi }_i={\lambda }_i\mathbf{M}{\varphi }_i$
请注意，Vallet & Levy 提出解决对称标准特征值问题的方法:
${\mathbf{M}}^{-1/2}\mathbf{L}{\mathbf{M}}^{-1/2}{\varphi }_i={\lambda }_i{\varphi }_i$
Libigl 实现了一个广义特征问题求解器，因此可以避免这种不必要的对称化。
通常上面的总和被截断为第一个$k$特征向量。如果选择了低频模式，即那些对应于小的$\lambda$值的模式，那么这种截断有效地将$\mathbf{f}$正则化为平滑，即网格上缓慢变化的函数。
模态分析和模型子空间在实时变形中被频繁使用。
在示例 306中，计算离散 Laplace-Beltrami 算子的前 5 个特征向量，并在甲虫顶部以伪彩色显示。使用 igl::wigs（镜像 MATLAB eigs）计算特征向量。将 5 个特征向量放入 U 的列中，并将特征值放入 S 的条目中：

SparseMatrix<double> L,M;
igl::cotmatrix(V,F,L);
igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,M);
Eigen::MatrixXd U;
Eigen::VectorXd S;
igl::eigs(L,M,5,igl::EIGS_TYPE_SM,U,S);