Chapter 2: Matrices and Linear Algebra¶

import igl
import scipy as sp
import numpy as np
from meshplot import plot, subplot, interact

import os
root_folder = os.getcwd()

Laplace equation¶

A common linear system in geometry processing is the Laplace equation:

$∆z = 0$

subject to some boundary conditions, for example Dirichlet boundary conditions (fixed value):

$\left.z\right|_{\partial{S}} = z_{bc}$

In the discrete setting, the linear system can be written as:

$\mathbf{L} \mathbf{z} = \mathbf{0}$

where $\mathbf{L}$ is the $n \times n$ discrete Laplacian and $\mathbf{z}$ is a vector of per-vertex values. Most of $\mathbf{z}$ correspond to interior vertices and are unknown, but some of $\mathbf{z}$ represent values at boundary vertices. Their values are known so we may move their corresponding terms to the right-hand side.

Conceptually, this is very easy if we have sorted $\mathbf{z}$ so that interior vertices come first and then boundary vertices:

$\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} \mathbf{0}_{in}\\ \mathbf{z}_{bc}\end{array}\right)$

The bottom block of equations is no longer meaningful so we’ll only consider the top block:

$\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) = \mathbf{0}_{in}$

We can move the known values to the right-hand side:

$\mathbf{L}_{in,in} \mathbf{z}_{in} = - \mathbf{L}_{in,b} \mathbf{z}_{b}$

Finally we can solve this equation for the unknown values at interior vertices $\mathbf{z}_{in}$ .

However, our vertices will often not be sorted in this way. One option would be to sort V, then proceed as above and then unsort the solution Z to match V. However, this solution is not very general.

With array slicing no explicit sort is needed. Instead we can slice-out submatrix blocks ( $\mathbf{L}_{in,in}$ , $\mathbf{L}_{in,b}$ , etc.) and follow the linear algebra above directly. Then we can slice the solution into the rows of Z corresponding to the interior vertices.

from scipy.sparse.linalg import spsolve

v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "camelhead.off"))

## Find boundary vertices
e = igl.boundary_facets(f)
v_b = np.unique(e)

## List of all vertex indices
v_all = np.arange(v.shape[0])

## List of interior indices
v_in = np.setdiff1d(v_all, v_b)

## Construct and slice up Laplacian
l = igl.cotmatrix(v, f)
l_ii = l[v_in, :]
l_ii = l_ii[:, v_in]

l_ib = l[v_in, :]
l_ib = l_ib[:, v_b]

## Dirichlet boundary conditions from z-coordinate
z = v[:, 2]
bc = z[v_b]

## Solve PDE
z_in = spsolve(-l_ii, l_ib.dot(bc))

plot(v, f, z)

Quadratic energy minimization¶

The same Laplace equation may be equivalently derived by minimizing Dirichlet energy subject to the same boundary conditions:

$\mathop{\text{minimize }}_z \frac{1}{2}\int\limits_S \|\nabla z\|^2 dA$

On our discrete mesh, recall that this becomes

$\mathop{\text{minimize }}_\mathbf{z} \frac{1}{2}\mathbf{z}^T \mathbf{G}^T \mathbf{D} \mathbf{G} \mathbf{z} \rightarrow \mathop{\text{minimize }}_\mathbf{z} \mathbf{z}^T \mathbf{L} \mathbf{z}$

The general problem of minimizing some energy over a mesh subject to fixed value boundary conditions is so wide spread that libigl has a dedicated api for solving such systems.

Let us consider a general quadratic minimization problem subject to different common constraints:

$\mathop{\text{minimize }}_\mathbf{z} \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} + \mathbf{z}^T \mathbf{B} + \text{constant},$

subject to

$\mathbf{z}_b = \mathbf{z}_{bc} \text{ and } \mathbf{A}_{eq} \mathbf{z} = \mathbf{B}_{eq},$

where

$\mathbf{Q}$ is a (usually sparse) $n \times n$ positive semi-definite matrix of quadratic coefficients (Hessian),
$\mathbf{B}$ is a $n \times 1$ vector of linear coefficients,
$\mathbf{z}_b$ is a $|b| \times 1$ portion of $\mathbf{z}$ corresponding to boundary or fixed vertices,
$\mathbf{z}_{bc}$ is a $|b| \times 1$ vector of known values corresponding to $\mathbf{z}_b$ ,
$\mathbf{A}_{eq}$ is a (usually sparse) $m \times n$ matrix of linear equality constraint coefficients (one row per constraint), and
$\mathbf{B}_{eq}$ is a $m \times 1$ vector of linear equality constraint right-hand side values.

This specification is overly general as we could write $\mathbf{z}_b = \mathbf{z}_{bc}$ as rows of $\mathbf{A}_{eq} \mathbf{z} = \mathbf{B}_{eq}$ , but these fixed value constraints appear so often that they merit a dedicated function: min_quad_with_fixed.

Linear equality constraints¶

We saw above that min_quad_with_fixed in libigl provides a compact way to solve general quadratic programs. Let’s consider another example, this time with active linear equality constraints. Specifically let’s solve the bi-Laplace equation or equivalently minimize the Laplace energy:

$\Delta^2 z = 0 \leftrightarrow \mathop{\text{minimize }}\limits_z \frac{1}{2} \int\limits_S (\Delta z)^2 dA$

subject to fixed value constraints and a linear equality constraint:

$z_{a} = 1, z_{b} = -1$ and $z_{c} = z_{d}$ .

Notice that we can rewrite the last constraint in the familiar form from above:

$z_{c} - z_{d} = 0.$

Now we can assembly Aeq as a $1 \times n$ sparse matrix with a coefficient $1$ in the column corresponding to vertex $c$ and a $-1$ at $d$ . The right-hand side Beq is simply zero.

Internally, min_quad_with_fixed solves using the Lagrange Multiplier method. This method adds additional variables for each linear constraint (in general a $m \times 1$ vector of variables $\lambda$ ) and then solves the saddle problem:

$\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\, \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} + \mathbf{z}^T \mathbf{B} + \text{constant} + \lambda^T\left(\mathbf{A}_{eq} \mathbf{z} - \mathbf{B}_{eq}\right)$

This can be rewritten in a more familiar form by stacking $\mathbf{z}$ and $\lambda$ into one $(m+n) \times 1$ vector of unknowns:

$\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\, \frac{1}{2} \left( \mathbf{z}^T \lambda^T \right) \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) + \left( \mathbf{z}^T \lambda^T \right) \left( \begin{array}{c} \mathbf{B}\\ -\mathbf{B}_{eq} \end{array} \right) + \text{constant}$

Differentiating with respect to $\left( \mathbf{z}^T \lambda^T \right)$ reveals a linear system and we can solve for $\mathbf{z}$ and $\lambda$ . The only difference from the straight quadratic minimization system, is that this saddle problem system will not be positive definite. Thus, we must use a different factorization technique (LDLT rather than LLT): libigl’s min_quad_with_fixed automatically chooses the correct solver in the presence of linear equality constraints.

The following example first solves with just fixed value constraints (left: 1 and -1 on the left hand and foot respectively), then solves with an additional linear equality constraint (right: points on right hand and foot constrained to be equal).

v, f = igl.read_triangle_mesh(os.path.join(root_folder, "data", "cheburashka.off"))

## Two fixed points: Left hand, left foot should have values 1 and -1
b = np.array([4331, 5957])
bc = np.array([1., -1.])
B = np.zeros((v.shape[0], 1))

## Construct Laplacian and mass matrix
L = igl.cotmatrix(v, f)
M = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)
Minv = sp.sparse.diags(1 / M.diagonal())

## Bi-Laplacian
Q = L @ (Minv @ L)

## Solve with only equality constraints
Aeq = sp.sparse.csc_matrix((0, 0))
Beq = np.array([])
_, z1 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)

## Solve with equality and linear constraints
Aeq = sp.sparse.csc_matrix((1, v.shape[0]))
Aeq[0,6074] = 1
Aeq[0, 6523] = -1
Beq = np.array([0.])
_, z2 = igl.min_quad_with_fixed(Q, B, b, bc, Aeq, Beq, True)

## Normalize colors to same range
min_z = min(np.min(z1), np.min(z2))
max_z = max(np.max(z1), np.max(z2))
z = [(z1 - min_z) / (max_z - min_z), (z2 - min_z) / (max_z - min_z)]

## Plot the functions
p = subplot(v, f, z[0], shading={"wireframe":False}, s=[1, 2, 0])
subplot(v, f, z[1], shading={"wireframe":False}, s=[1, 2, 1], data=p)

# @interact(function=[('z0', 0), ('z1', 1)])
# def sf(function):
#     p.update_object(colors=z[function])

/home/nico/git/libigl-python-bindings/venv/lib64/python3.8/site-packages/scipy/sparse/_index.py:103: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.
  self._set_intXint(row, col, x.flat[0])