介绍

在进行软体模拟时，如果使用牛顿法计算最优的下降方向，需要计算能量密度函数\(\Psi\)关于位置\(\vec {\mathbf{x}}\) 的Hessian矩阵，即\(\frac{\partial^2 \Psi}{\partial \vec {\mathbf{x}}^2}\)。其中\(\vec {\mathbf{x}}\)是一个四面体的四个顶点的位置。因为：

\[ \begin{aligned} \frac{\partial^2 \Psi}{\partial \vec {\mathbf{x}}^2} &= \frac{\partial}{\partial \vec {\mathbf{x}}} \left(\frac{\partial \Psi}{\partial \vec {\mathbf{x}}}\right) \\ &= \frac{\partial}{\partial \vec {\mathbf{x}}} \left(\frac{\partial \Psi}{\partial F}:\frac{\partial F}{\partial \vec {\mathbf{x}}}\right) \\ &= \frac{\partial}{\partial \vec {\mathbf{x}}} \left(\frac{\partial \Psi}{\partial F}\right):\frac{\partial F}{\partial \vec {\mathbf{x}}} + \frac{\partial \Psi}{\partial F}:\frac{\partial^2 F}{\partial \vec {\mathbf{x}}^2}\\ &=\frac{\partial}{\partial \vec {\mathbf{x}}} \left(P\right):\frac{\partial F}{\partial \vec {\mathbf{x}}}\\ &=(\frac{\partial F}{\partial \vec {\mathbf{x}}}: \frac{\partial P}{\partial F}):\frac{\partial F}{\partial \vec {\mathbf{x}}} \end{aligned} \]

证明

首先，已知

\[ \vec {\mathbf{x}} = \begin{bmatrix} \vec x_1 & \vec x_2 & \vec x_3 \end{bmatrix}\\ D_s = \begin{bmatrix} \vec x_1-\vec x_4 & \vec x_2-\vec x_4 & \vec x_3-\vec x_4 \end{bmatrix}\\ \begin{aligned} &\frac{\partial (D_s)_{kl}}{\partial \vec {\mathbf{x}_{ij}}} e_i\otimes e_j\otimes e_k\otimes e_l\\ &=\frac{\partial \vec {\mathbf{x}}_{kl}}{\partial \vec {\mathbf{x}}_{ij}} e_i\otimes e_j\otimes e_k\otimes e_l\\ &= \delta_{ik}\delta_{jl} e_i\otimes e_j\otimes e_k\otimes e_l \end{aligned} \] \(D_m^{-1}\)的分量表示为\(d_{mn}\)，\(P\)的分量表示为\(P_{rs}\)，则能量密度函数\(\Psi\)关于位置\(\vec {\mathbf{x}}\) 的梯度为：

Understading the Geometry of \(A\in \mathbb{R}^{m\times n}\)

\[ A\in \mathbb{R}^{m\times n} \\ x\in \mathbb{R}^n \mapsto Ax \in \mathbb{R}^m \]

Assumption: \(A \in \mathbb{R}^{n \times n}\) is symmetric, no repeated eigenvalues \(\lambda\)s.

Two basic Properties

Lemma: Every matrix \(A \in \mathbb{R}^{n \times n}\) has at least one (complex) eigenvalue.

很多时候需要把已有问题转化为特征值问题。

1. 有一个对称矩阵（有时候是对称半正定矩阵），我们想要找一个方便的基。
2. 具体的优化问题。

High-Level Idea

Why QR?

\[ \begin{aligned} cond A^T A &= ||A^T A|| \cdot ||(A^T A)^{-1}||\\ &\approx ||A^T|| \cdot ||A|| \cdot ||A^{-1}|| \cdot ||A^{-T}|| &= cond A^2 \end{aligned} \]

为了避免计算\(A^T A\)，我们可以使用QR分解。

介绍

Projective Dynamics是一种用于软体模拟的方法，算法分为local step和global step两个部分。其中local step可以对于每个四面体约束并行计算，global step只需要求解一个线性方程组，而他的矩阵非常特殊，是一个Gram矩阵，因此可以预先用Cholesky分解。本文主要介绍local step的实现的推导过程。用Corotated strain model的CUDA的实现作为例子。

Sensitivity and Conditioning

Gaussian elimination works in theory, but what about floating point precision?

How much can we trust \(x_0\) if \(0<||Ax_0 - b||_2 \ll 1\)(backwards error)?

Structured Linear Systems

Gaussian elimination and/or LU can solve all the example problems above. But these systems can have special properties that make them easier or stabler to solve.

Today's example: Positive definite, sparsity.