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

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

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

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:

subject to

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:

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:

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

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])
```