Getting to quadratic form

[oveddan]

I'm not sure really how to go from:

to:

If there was no A, it could be done like (Where Gt = G transposed): GtGs=Gtu

Are there any examples out there for how this could be done with an additional matrix such as A outside of the least squares term?

[jpanetta]

I think this should be straightforward if you understand where G^T G comes from.

Pretend for the moment that the energy doesn't have area weights: E_1(s) := sum_t ||g_t - u_t||^2 = sum_t ||G_t s - u_t||^2 where G_t is a 3x#V matrix mapping the per-vertex scalar values to the gradient on triangle t.

Expanding this energy, we see E_1(s) = sum_t s^T G_t^T G_t s - 2 s^T G_t^T u_t + ||u_t|| = s^T (sum_t G_t^T G_t) s - 2 s^T (sum_t G_t^T u_t) + const

Now, the "G" matrix I believe you refer to, which can be computed by igl::grad, holds the G_t for each triangle t (but not directly stacked: it chooses to place row i of G_t in G.row(i * #F + t)). You should be able to see that G^T G is equivalent to (sum_t G_t^T G_t): the matrix multiplication implements the sum over triangles (as well as over the rows of per-triangle G_t matrices).

The same can be done with the linear term as long as you flatten u_t into a column vector "u" using the same ordering as "G" (i.e. place row i of u_t in u.row(i * #F + t)): sum_t G_t^T u_t = G^T u

So the energy looks like: E_1(s) = s^T G^T G s - 2 s^T G^T u + const

Now, if we were to include the triangle areas, we'd have: E_A(s) = s^T (sum_t G_t^T A_t G_t) s - 2 s^T (sum_t G_t^T A_t u_t) + const, which you can interpret as scaling one of the G_t in each of the sums by A_t. This in turn means one of the G matrices in each term of "s^T G^T G s - 2 s^T G^T u" should have its rows scaled by the area of the triangles from which they originated. Finally, recall that you can scale the rows of a matrix by applying a diagonal matrix on the left--this should be all you need.

Implementation note: diagonal matrices can be implemented efficiently in Eigen by constructing a vector holding the diagonal entries and calling "asDiagonal()."

[oveddan]

Thank you, this makes more sense now.