Minimizing quadratic energies with constant constraints

Alec Jacobson

June 21, 2011

weblog/

There are a many different ways to minimize a quadratic energy with constant constraints. I wrote out the exact steps for three such ways:

Separate the energy into knowns and unknowns, then solve for the unknowns by setting gradient with respect to those unknowns to zero.
Add additional term to energy that punishes not satisfying the known values of your constant constraints in the least squares sense, the solve for all variables by setting gradient with respect to all variables to zero.
Add additional variables called lagrange multipliers to energy so that the problem is no longer to find a minimum but a stationary point in the new energy. Such a stationary point guarantees that our constraints are met. Find the stationary point by the usual method of setting the gradient to zero, this time the gradient is with respect to all original variables and these new variables.

The pretty print type-up with MATLAB pseudocode. A plain text, unicode version:

Quadratic energy
E = ZT * Q * Z + ZT * L + C

We want to fix some of the values of Z, we can rearrange the Z's that are still
free and the ones we want to fix so that Z = [X;Y], where X are the free values
and Y are the fixed values;

Split energy into knowns and unknowns
Z = [X;Y]
or 
Z(unknown) = X and Z(known) = Y
E = [X;Y]T * Q * [X;Y] + [X;Y]T * L + C
E = [X;0]T * Q * [X;Y] + [0;Y]T * Q * [X;Y] + [X;0] * L + [X;Y]T * L + C
E = [X;0]T * Q * [X;0] + 
    [X;0]T * Q * [0;Y] + 
    [0;Y]T * Q * [X;0] + 
    [0;Y]T * Q * [0;Y] + 
    [X;0]T * L + 
    [0;Y]T * L + 
    C
E = [X;0]T * [Q(unknown,unknown) 0; 0 0] * [X;0] + 
    [X;0]T * [0 Q(unknown,known); 0 0] * [0;Y] + 
    [0;Y]T * [0 0; Q(known,unknown) 0] * [X;0] + 
    [0;Y]T * [0 0; 0 Q(known,known)] * [0;Y] + 
    [X;0]T * L + 
    [0;Y]T * L + 
    C
E = XT * Q(unknown,unknown) * X +
  XT * Q(unknown,known) * Y +
  YT * Q(known,unknown) * X +
  YT * Q(known,known) * Y + 
  XT * L(unknown) +
  YT * L(known) +
  C

E = XT * Q(unknown,unknown) * X +
  XT * Q(unknown,known) * Y +
  (YT * Q(known,unknown) * X)T +
  YT * Q(known,known) * Y + 
  XT * L(unknown) +
  YT * L(known) +
  C

E = XT * Q(unknown,unknown) * X +
  XT * Q(unknown,known) * Y +
  XT * Q(known,unknown)T * Y +
  YT * Q(known,known) * Y + 
  XT * L(unknown) +
  YT * L(known) +
  C

E = XT * NQ * X + XT * NL + NC
NQ = Q(unknown,unknown)
NL = Q(unknown,known) * Y + Q(known,unknown)T * Y + L(unknown)
NC = YT * Q(known,known) * Y + YT * L(known) + C

∂E/∂XT = 2 * NQ * X + NL

Solve for X with:
X = (- 2 * NQ)^-1 * NL

Enforce fixed values via soft contraints with high weights
E = ZT * Q * Z + ZT * L + C 
Add new energy term punishing being far from fixed variables
NE = ZT * Q * Z + ZT * L + C + w * (Z(known) - Y)T * I * (Z(known) - Y)
where w is a large number, e.g. 10000
NE = ZT * Q * Z + ZT * L + C + w * (Z - [0;Y])T * [0 0;0 I(known,known)] * (Z - [0;Y])
where W is a diagonal matrix with Wii = 0 if i∈unknown and Wii = wii if j∈known
NE = ZT * Q * Z + ZT * L + C + (Z - [0;Y])T * W * (Z - [0;Y])
NE = ZT * Q * Z + ZT * L + C + ZT * W * Z - 2 * ZT * W * [0;Y] + [0;Y]T * W * [0;Y]
NE = ZT * NQ * Z + ZT * NL + NC
NQ = Q + W
NL = L - 2W * [0;Y]
NC = C + [0;Y]T * W * [0;Y]

Differentiate with respect to ZT
∂E/∂ZT = 2 * NQ * Z + NL

Solve for Z with:
Z = -0.5 * NQ^-1 * NL
or re-expanded
Z = -0.5 * (Q + W)^-1 * (L - 2 * W * [0;Y])

Discard known parts of Z to find X
X = Z(unknown)
But be careful to look at Z(known) - Y to be sure your fixed values are being
met, if not then w is too low. If the X values look like garbage then perhaps
you're getting numerical error because w is too high.

Lagrange multipliers
We want to minimize
E = ZT * Q * Z + ZT * L + C 
subject to the constant constraints:
Z(known) = Y
Find stationary point of new energy:
NE = E + ∑λi *(Zi - Yi)
        i∈known
NE = E + [Z;λ]T * QC * [Z;λ] - [Z;λ]T * [0;Y]
where QC = [0 0 0;0 0 I(known,known);0 I(known,known) 0]
NE = [Z;λ]T * NQ * [Z;λ] + [Z;λ]T * NL + C 
NQ = [Q 0;0 0] + QC =  [         Q        0             ]
                       [                  I(known,known)]
                       [ 0 I(known,known) 0             ]
NL = [L;0] + [0;Y];
Differentiate with respect to all variables, including lagrange multipliers
∂E/∂[Z;λ]T  = 2 * NQ * Z + NL
Solve with
[Z;λ] = -0.5 * NQ^-1 * NL

Discard fixed values and langrange multipliers.
X = Z(known)
The value of λi is the "force" of the constraint or how hard we had to pull the
energery minimum to meet the constraints. See
http://www.slimy.com/~steuard/teaching/tutorials/Lagrange.html#Meaning