Questions about Preconditioning Conjugate Gradient method in Baraff & Witkin 98?

Question

I encountered some problems when implementing the cloth simulation algorithm from Baraff & Witkin 98's Large Steps in Cloth Simulation.

---

Baraff & Witkin 98

Consider the cloth as a particle system in 3D, which consists of N particles.

1. Equation of motion

In each time step, the implicit integration of particle system is:

$$ \begin{pmatrix} \Delta \mathbf x \\ \Delta \mathbf v \\ \end{pmatrix} = h\begin{pmatrix} \mathbf v_0 + \Delta \mathbf v \\ \mathbf M^{-1} \mathbf f(\mathbf x_0 + \Delta \mathbf x, \mathbf v_0 + \Delta \mathbf v) \\ \end{pmatrix} $$

where $h$ is the time step (scalar)

2. The partial derivative equation

Rewrite the aforementioned equality of motion into $\mathbf A \Delta \mathbf v = \mathbf b$ form, ( Take the 1st order Taylor expansion of $\mathbf f(\mathbf x_0 + \Delta \mathbf x, \mathbf v_0 + \Delta \mathbf v)$, and introduce constraint matrix $\mathbf S$, which make mass matrix from $\mathbf M^{-1}$ to $\mathbf M^{-1}\mathbf S$), we have:

$$ (\mathbf I - h\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf v} - h^2\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf x})\Delta \mathbf v = h\mathbf M^{-1} \mathbf S(\mathbf f_0 + h\frac{\partial \mathbf f}{\partial \mathbf x} \mathbf v_0) + \mathbf z $$

Where all of the matrices are treated as $N\times N $ block matrices, eack block sizes $3\times 3 $; vectors as $N$ block vector, each sized $3$ (with N number of 3D forces/positions/velocities).

• $\mathbf I$ is identical matrix

• $ \mathbf M = diag(\mathbf M_1, \mathbf M_2, ..., \mathbf M_N) $, $\mathbf M_i = diag(m_i, m_i, m_i)$, mass of each particle

• $ \mathbf S = diag(\mathbf S_1, \mathbf S_2, ..., \mathbf S_N) $, $ \mathbf S_i = \begin{cases} \mathbf I, & \text{$ndof(i) = 3$} \\ \mathbf I - \mathbf p_i \mathbf p_i^T, & \text{$ndof(i) = 2$} \\ \mathbf I - \mathbf p_i \mathbf p_i^T - \mathbf q_i \mathbf q_i^T, & \text{$ndof(i) = 1$} \\ \mathbf 0, & \text{$ndof(i) = 0$} \end{cases}$, $\mathbf p$ and $\mathbf q$ are two oithogonal constraint direction, $\mathbf S$ is constraint matrix (Baraff98 chapter 5.1)

• $\frac{\partial \mathbf f}{\partial \mathbf v}$ and $\frac{\partial \mathbf f}{\partial \mathbf x}$ are force derivatives, which are symmetrical matrices

.

3. solve the PDE

The paper solves the aforementioned PDE by Modified Preconditioning Conjugate Gradient Method (Baraff98 chapter 5.2) for $\Delta \mathbf v$, then update the velocity $\mathbf v$ and position $\mathbf x$ of each particle.

---

My Questions

Solve $\mathbf A\mathbf x=\mathbf b$ with PCG method requires the matrix $\mathbf A$ to be symmetrical and positive-definite, where according to the aforementioned PDE, $\mathbf A = (\mathbf I - h\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf v} - h^2\mathbf M^{-1} \mathbf S\frac{\partial \mathbf f}{\partial \mathbf x})$.

1. symmentrical

In block view: treat $\mathbf A$ as an $N\times N$ block matrix, $\mathbf M^{-1}$ and $\mathbf S$ are both block diagonal matrices, makes $\mathbf A$ a block symmetrical matrix. While in normal view: $\mathbf A$ as $3N\times 3N$ matrix, $\mathbf S$ is symmetrical, but product of $\mathbf M^{-1} \mathbf S$ and $\frac{\partial \mathbf f}{\partial \mathbf v}$ (or $\frac{\partial \mathbf f}{\partial \mathbf x}$) makes $\mathbf A$ not symmetrical.

Question 1: should I treat $\mathbf A\mathbf x=\mathbf b$ as a block matrix rather than normal matrix? If not, how to make $\mathbf A$ symmetrical?

I then modified the PDE as $\mathbf A^T\mathbf A\mathbf x=\mathbf A^T\mathbf b$, $\mathbf A^T\mathbf A$ is symmetrical.

Question 2: How to make the fliter procedure (Baraff98 chapter 5.3) compatible with $\mathbf A^T\mathbf A\mathbf x=\mathbf A^T\mathbf b$?

2. positive-definite

The blocks $\mathbf S_i$ in $\mathbf S$ may become zero block when $ndof(i) = 0$, which makes $\mathbf A^T\mathbf A$ a positive-semidefinite matrix. To apply PCG method, according to Baraff98 chapter 5.3:

The CG method (technically, the preconditioned CG method) takes a symmetric positive semi-definite matrix $\mathbf A$, a symmetric positive definite preconditioning matrix $\mathbf P$ of the same dimension as $\mathbf A$ ,a vector $\mathbf b$ and iteratively solves $\mathbf A \Delta \mathbf v = \mathbf b$.

Question 3: How to find the symmetrical and positive-definite preconditioning matrix $\mathbf P$, with $\mathbf A^T\mathbf A$ symmetrical and semi PD?

---

Some docs I also refered:

• On the modified conjugate gradient method in cloth simulation
• freecloth
• Cloth Simulation slides.

Answer 1

Baraff and Witkin propose to incorporate constraints by modifying the linear system by the constraints matrices $\mathbf S_i$. As they state in the beginning of Section 5.2, the resulting system is not symmetric anymore. Therefore, the modified linear system must be computed with a linear solver that can treat non-symmetric systems (as - if I got you correctly - you proposed with the normal matrix $\mathbf A^T \mathbf A$). However, Baraff and Witkin introduce the modified conjugate gradient method, which is a variant of the original conjugate gradient method, that computes a solution of the initial linear system with considering the constraints at the same time, without the need to take care of the non-symmetry. In their paper they do not solve the normal equations ($\mathbf A^T \mathbf A$), but rather the original system $\mathbf A \Delta \mathbf v = \mathbf b$, so your questions ask something which is beyond of the paper.