Lecture 12 Deformable Simulation 02|The Implicit Integration Methods

# The Implicit Euler Integration

# Time Step

• In simulation: The time difference between the adjacent ticks on the temporal axis for the simulation: $h_{\text{sim}}$
  • v[i] += h * acc
  • x[i] += h * v[i]
• In display: The time difference between two images displayed on the screen: $h_{\text{disp}} = 1/60\ \mathrm{s}$ for 60 Hz applications

Sub-(time)-stepping: $n_{\text{sub}}$: $h_{\text{sim}} = \frac{h_{\text{disp}}} {n_{\text{sub}}}$

The smaller $n_{\text{sub}}$, the larger $h_{\text{sim}}$ (less accurate)

• For explicit integrations: too long time step may lead to explosion (as $h$ increases, diverge);
• For implicit: usually converge for all $h$ => artifact: slower converge

# Numerical Recipes for Implicit Integrations

# The Baraff and Witkin Style Solution

One iter of Newton’s Method, referred as semi-implicit

• $x_{n+1} = x_n + hv_{n+1}$ & $v_{n+1} = v_n + hM^{-1} f(x_{n+1})$ then $\Rightarrow x_{n+1} = x_n + hv_n + h^2 M^{-1}f(x_{n+1})$

• Goal: Solving $x_{n+1} = x_n + hv_n + h^2 M^{-1}f(x_{n+1})$

• Assumption: $x_{n+1}$ is not too far away from $x_n$

• Algorithm: Non-linear -> Linear ($\delta x$ - the deviation between the two timesteps, which is very small)

  • Let $\delta x = x_{n+1} - x_n$, then f(x_{n+1}) \approx f(x_n) + \grad_{x}f(x_n)\delta x (1st order Taylor, neglect all the minor errors)

  • Substitute this approx:

    • $x_{n}+\delta x=x_{n}+h v_{n}+h^{2} M^{-1}\left(f\left(x_{n}\right)+\nabla_{x} f\left(x_{n}\right) \delta x\right)$
    • $\left(M-h^{2} \nabla_{x} f\left(x_{n}\right)\right) \delta x=h M v_{n}+h^{2} f\left(x_{n}\right)$ (Matrix op: (2n x 2n) x (2n x 1) = (2n x 1), sim to Ax = b)
    • $x_{n+1} = x_n + \delta x$, $v_{n+1} = \delta x/h$

    The solution is actually the location of $x_{n+1}$ while in implicit method the sol will be much more closed to the exact root (error is the red line)

    

# Reformulating the Implicit Euler Problem

To reduce the red part: add another step:

• Integrating the nonlinear root finding problem over $x$ gives:

  x_{n+1}=\operatorname{argmin}_{x}\left(\frac{1}{2}\left\|x-\left(x_{n}+h v_{n}\right)\right\|_{M}^{2}+h^{2} E(x)\right), \text { given } f(x)=\nabla_{x} E(x)

  Note: (Matrix Norm) $\| x\| _A = \sqrt{x^{\mathrm{T}}Ax}$; (Vector Derivative) \grad_x (x^{\mathrm{T}}Ax) = (A +A^{\mathrm{T}})x

# Minimization / Non-linear Root-finding:

• Let $g(x) =\frac{1}{2}\left\|x-\left(x_{n}+h v_{n}\right)\right\|_{M}^{2}+h^{2} E(x)$ (Mass matrix is usually diagonalised, so the transpose is itself)

• Then $\nabla_{x} g\left(x_{n+1}\right)=M\left(x_{n+1}-\left(x_{n}+h v_{n}\right)\right)-h^{2} f\left(x_{n+1}\right)$ (Energy of the der of pos => The fastest energy increasing dir, the negative val is force)

• For nonsingular $M$: ($g$ at an extreme => the derivatifve of $g$ at 0, the equation above equals to 0)

  $\nabla_{x} g\left(x_{n+1}\right)=0 \leftrightarrow x_{n+1}=\left(x_{n}+h v_{n}\right)+h^{2} M^{-1} f\left(x_{n+1}\right)$ (The root of implicit Euler)

# Convex Minimization of $\min g(x)$

The general descent method:

def minimize_g():
    x = x_0
    while grad_g(x).norm() > EPSILON:
        Determine a descent dir: dx
        Line search: choose a stepsize: t > 0
        Update: x = x + t * dx

• Determine a descent dir: $\mathrm{d}x $

  • Opt 1: Gradient Descent: \mathrm{d}x = -\grad_x g(x) (use the intersect point with the x-axis)

    

  • Opt 2: Newton’s Method: $\mathrm{d} x=-\left(\nabla_{x}^{2} g(x)\right)^{-1} \nabla_{x} g(x)$ (apply the second order derivative as the curvature -> second order -> use its valley (actually also decide the step)) -> usually when a func can solve 2nd order

    

• Find the step size:

  • Line search: choose a step size $t>0$ given a $\mathrm{d}x$

  • Backtracking:

    def line_search(): 
        t = t_0
        while g(x) + alpha*t*grad_g(x).dot(dx) < g(x+t*dx): 
            # this alpha indicates a line between the grad and the line w/ alpha = 0
            t = t * beta
            
        return t
    

    

    $\alpha \in (0, 0.5)$, $\beta \in (0,1)$ (common choice: $\alpha = 0.03 $ & $\beta = 0.5$)

Problem

But most deformable bodies have non-convex energies: unsteady equilibrium.

# Steps

$g(x)$ is contains 2 terms: dynamics and elastic

$g(x)=\frac{1}{2}\left\|x-\left(x_{n}+h v_{n}\right)\right\|_{M}^{2}+h^{2} E(x)$

1. Init guess: $x = x_n$ or $x = x_n + hv_n$
2. While loop: while not converge:
   • Descent dir: \mathrm{d}x = -\grad_x g(x) or \mathrm{d}x = -(\grad_x^2 g(x))^{-1} \grad_x g(x)
     • Gradient: \grad_{x} g(x)=M\left(x-\left(x_{n}+h v_{n}\right)\right)+h^{2} \grad_{x} E(x) (linear + non-linear) (Model dep)
     • Hessian (Matrix): \grad_x^2 g(x= M+ h^2\grad^2_x E(x) (Model dep)
   • Line search: det the stepsize $t$
   • Update: $x = x + t\cdot \mathrm{d}x$

# Definiteness-fixed Newton’s Method

For general cases:

in Step 2 (while), After compute gradient dir and Hessian, add another substep:

• Fix Hessian to positive definite: $\tilde{H} = \mathrm{fix}(H)$

# Definiteness-fix

$\tilde{H} = \mathrm{fix}(H)$

• Opt 1: Global Regularization

  • Init: $\tilde{H} = H$, flag = False, reg = 1
  • while not flag:
    • flag, L = factorize(~H) Try to factorize $\tilde{H} = LL^{\mathrm{T}}$ (if success -> OK; else: add identity to let it tends to indentity (definite))
    • $\tilde{H} = H + \mathrm{reg} * I$, reg = reg * 10

• Opt 2: Local Regularization (in mass-spring system)

  $$ \nabla_{x}^{2} E(x)=\left[\begin{array}{cccc} K_{1} & -K_{1} & & \\ -K_{1} & K_{1}+K_{2}+K_{3} & -K_{3} & -K_{2} \\ & -K_{3} & K_{3} & \\ & -K_{2} & & K_{2} \end{array}\right] $$

  $\left[\begin{array}{cc} K_{1} & -K_{1} \\ -K_{1} & K_{1} \end{array}\right]\left[\begin{array}{cc} K_{2} & -K_{2} \\ -K_{2} & K_{2} \end{array}\right]\left[\begin{array}{cc} K_{3} & -K_{3} \\ -K_{3} & K_{3} \end{array}\right]$

  • $\nabla_{x}^{2} g(x)=M+h^{2} \nabla_{x}^{2} E(x)=M+h^{2} \sum_{j=1}^{m} \nabla_{x}^{2} E_{j}(x)$

    • $K_{1} \geqslant 0 \Rightarrow\left[\begin{array}{cc} K_{1} & -K_{1} \\ -K_{1} & K_{1} \end{array}\right]=K_{1} \otimes\left[\begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array}\right] \geqslant 0$ (a positive semi-definite matrix Otimes another .. matrix => also positive semi-definite)
    • $\Rightarrow \nabla_{x}^{2} E(x)=\nabla_{x}^{2} E_{1}(x)+\nabla_{x}^{2} E_{2}(x)+\nabla_{x}^{2} E_{3}(x) \geqslant 0$
    • $ \Rightarrow \nabla_{x}^{2} g(x)=M+h^{2} \nabla_{x}^{2} E(x)>0$

    Has a sufficient condition: $K_1 \ge 0,\ K_2 \ge0, \ K_3\ge0$

  • $K_{1}=k_{1}\left(I-\frac{l_{1}}{\left\|x_{1}-x_{2}\right\|}\left(I-\frac{\left(x_{1}-x_{2}\right)\left(x_{1}-x_{2}\right)^{T}}{\left\|x_{1}-x_{2}\right\|^{2}}\right)\right) \in \mathbb{R}^{2 \times 2}$

    • $K_1 = Q \Lambda Q^{T}$ (<- Eigen value decomposition)
    • $\tilde \Lambda = \max(0,\Lambda)$ (<- Clamp the negative eigen values)
    • $\tilde{K_1} = Q\tilde{\Lambda} Q^{\mathrm{T}}$ (<- Construct the p.s.d. projection)

After this definiteness-fix: (Newton’s)

# Linear Solvers

# Linear Solvers for $Ax = b$

• Direct solvers

  • Inversion: $x = A^{-1} b$ (requires higher complexity for sparse and higher order matrices)

  • Factorization (usually better especially for sparse matrices):

    $ A=\left{\begin{array}{cl} L U & \text {, if } A \text { is a square matrix } \ L D L^{\mathrm{T}} & \text {, if } A=A^{T} \ L ^{\mathrm{T}} & \text {, if } A=A^{T} \text { and } A>0 \end{array}\right.$

• Interative solvers

  • Stationary iterative linear solvers: Jacobi / Gauss-Seidel / SOR / Multigrid
  • Krylov subspace methods: Conjugate Gradient (CG) / biCG / CR / MinRes / GMRes

# Factorization

For $A = L L^{\mathrm{T}}$ (Assume already p.s.d. proceeded)

• Solve $Ax = b$ is equiv to $LL^{\mathrm{T}}x = b$

• First solve $Ly = b$ (Forward substitution)

• Then solve $L^{\mathrm{T}} x = y$ (Backward substitution)

  For sparse matrices -> complexity not high; But not parallized (CPU backend only)

Current APIs:

• [SparseMatrixBuilder] = ti.linalg.SparseMatrixBuilder()
• [SparseMatrix] = [SparseMatrixBuilder].build()
• [SparseMatrixSolver] = ti.linalg.SparseSolver(solver_type, ordering)
• [NumpyArray] = [SparseMastrixSolver].solve([Field])

# Conjugate Gradient (CG) Method

Properties:

• Works for any symmetric positive definite matrix $A = A^{\mathrm{T}},\ A >0$
• Guarantees to converge in $n$ iter for A\in\R ^{n\times n}
• Works amazingly good if the condition number $\kappa = \lambda_{\max} / \lambda_{\min}$ (eigenvalues) of $A$ is small (when $\kappa = 1$, identity matrix)
• Short codes

Demo: (Python scope - control flow)

def conjugate_gradient(A, b, x):
	i = 0
	r = b – A @ x
	d = r
	delta_new = r.dot(r)
	delta_0 = delta_new
	while i < i_max and delta_new/delta_0 > epsilon**2:
		q = A @ d	# Matrix-Vector multiplication -> expensive
		alpha = delta_new / d.dot(q)	# dot prod -> expensive
		x = x + alpha*d
		r = b - A @ x # r = r - alpha * q
		delta_old = delta_new
		delta_new = r.dot(r)
		beta = delta_new / delta_old
		d = r + beta * d
		i = i + 1
	return x

# Accelerate the CG Method

• Reduce the time of sparse-matrix-vector multiplication of q = A @ d: (Piece-wise multiply)

  $A =M +h^{2} \sum_{j=1}^{m} \nabla_{x}^{2} E_{j}(x) \quad \Rightarrow A d=M d+h^{2} \sum_{j=1}^{m} \nabla_{x}^{2} E_{j}(x) d$

  Use @ti.kernel -> computing

  (Taichi enables thread local storage automatically for this reduction problem) (Taichi TLS / CUDA Reduction Guide)

• Reduce the condition number of $A$:

  The error (Upper bound)

  $\left\|e_{i}\right\|_{A} \leq 2\left(\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}\right)^{i}\left\|e_{0}\right\|_{A}$

• Instead of solving $Ax = b$ => solve for $M^{-1} Ax = M^{-1} b$ ($M\sim A$ as much as possible)

  • Jacobi: $M = \mathrm{diag}(A)$
  • Incomplete Cholesky: $M = \tilde{L}\tilde{L}^{\mathrm{T}}$
  • Multigri