<p>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.</p>

CC BY-SA 3.0
2016-02-16T01:14:47.317Z
Follow your answer, I found out that the Modified Conjugate Gradient method is proposed to deal with the non-symmetric semi-positive-definite linear equation(section 5.2). Thanks a lot.
stanleyerror
CC BY-SA 3.0
2016-02-15T21:30:17.253Z
2016-02-15T21:30:17.253Z
dweber
<p>I encountered some problems when implementing the cloth simulation algorithm from Baraff &amp; Witkin 98's <a href="http://www.cs.cmu.edu/~baraff/papers/sig98.pdf" rel="nofollow">Large Steps in Cloth Simulation</a>.</p>

<hr>

<p><strong>Baraff &amp; Witkin 98</strong></p>

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

<ol>
<li>Equation of motion</li>
</ol>

<p>In each time step, the implicit integration of particle system is:</p>

<p>$$
\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}
$$</p>

<p>where $h$ is the time step (scalar)</p>

<ol start="2">
<li>The partial derivative equation</li>
</ol>

<p>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: </p>

<p>$$
(\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
$$ </p>

<p>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).</p>

<ul>
<li><p>$\mathbf I$ is identical matrix</p></li>
<li><p>$ \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</p></li>
<li><p>$ \mathbf S = diag(\mathbf S_1, \mathbf S_2, ..., \mathbf S_N) $, 
$  \mathbf S_i = \begin{cases}  \mathbf I,  &amp; \text{$ndof(i) = 3$} \\ \mathbf I - \mathbf p_i \mathbf p_i^T,  &amp; \text{$ndof(i) = 2$} \\ \mathbf I - \mathbf p_i \mathbf p_i^T - \mathbf q_i \mathbf q_i^T,  &amp; \text{$ndof(i) = 1$} \\ \mathbf 0,  &amp; \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)</p></li>
<li><p>$\frac{\partial \mathbf f}{\partial \mathbf v}$ and $\frac{\partial \mathbf f}{\partial \mathbf x}$ are force derivatives, which are symmetrical matrices</p></li>
</ul>

<p>.</p>

<ol start="3">
<li>solve the PDE</li>
</ol>

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

<hr>

<p><strong>My Questions</strong></p>

<p>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})$.</p>

<ol>
<li>symmentrical</li>
</ol>

<p>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. </p>

<p><strong>Question 1</strong>: 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?</p>

<p>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.</p>

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

<ol start="2">
<li>positive-definite</li>
</ol>

<p>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:</p>

<blockquote>
  <p>The CG method (technically, the preconditioned CG method) takes a symmetric <strong>positive semi-definite</strong> matrix $\mathbf A$, a <strong>symmetric positive definite</strong> 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$.</p>
</blockquote>

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

<hr>

<p>Some docs I also refered:</p>

<ul>
<li><a href="http://www.cs.ubc.ca/~ascher/papers/ab.pdf" rel="nofollow">On the modified conjugate gradient method in cloth simulation</a></li>
<li><a href="http://davidpritchard.org/freecloth/docs/report-single/" rel="nofollow">freecloth</a></li>
<li><a href="http://cg.informatik.uni-freiburg.de/course_notes/sim_03_cloth1.pdf" rel="nofollow">Cloth Simulation slides</a>.</li>
</ul>

<real-time><mesh><physically-based><animation><simulation>
Questions about Preconditioning Conjugate Gradient method in Baraff & Witkin 98?
stanleyerror
# Questions about Preconditioning Conjugate Gradient method in Baraff & Witkin 98?

(Asked by: username_0 on 2016-02-02)

I encountered some problems when implementing the cloth simulation algorithm from Baraff & Witkin 98's [Large Steps in Cloth Simulation](http://www.cs.cmu.edu/~baraff/papers/sig98.pdf).

---

**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](http://www.cs.ubc.ca/~ascher/papers/ab.pdf)
* [freecloth](http://davidpritchard.org/freecloth/docs/report-single/)
* [Cloth Simulation slides](http://cg.informatik.uni-freiburg.de/course_notes/sim_03_cloth1.pdf).

**Tags:** real-time, mesh, physically-based, animation, simulation


---

## 1 answers

### Answer by username_1 on 2016-02-15. Score: 3 | Accepted answer

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.

> **Comments:**
> **username_0**: Follow your answer, I found out that the Modified Conjugate Gradient method is proposed to deal with the non-symmetric semi-positive-definite linear equation(section 5.2). Thanks a lot.  

---

