Skip to content

Chapter 2: Matrices and Linear Algebra

PyPI version buildwheels

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", ""))

## 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,

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},


  • \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", ""))

## 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/ SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil_matrix is more efficient.
  self._set_intXint(row, col, x.flat[0])