Ytukyg

I have an energy I'd like to minimize of the form:

$E(G) = |underbrace{X - Y G^T L G B}_{f(G)}|_F^2$

where $X,Y,B$ are dense matrices and $G,L$ are sparse matrices ($G^TLG$ is also sparse), $|M|_F^2$ computes the squared Frobenius norm of the matrix $M$.

I'd like to compute quantities such as $partial f/partial G$ (or $partial^2 E/partial G^2$) to use Gauss-Newton's (or Newton's) method, but I only care about changes to $G$ that maintain its sparsity pattern. That is, suppose $G$ is a function of the non-zero values collected in a vector $g$ via, say, $G = mathop{sparse}(i,j,g)$ where $i,j$ are lists of subscript indices corresponding to values in $g$. Then I'd really like to compute $partial f/partial g$.

Eventually I'm programming this and would like to avoid computing $partial f/partial G$ as a dense matrix/tensor.

Is there's a nice (reduced) expression for $partial f/partial g$?

If not, what's the best way to use automatic differentiation to conduct this optimization? I'm having trouble formulating this in the right way to use the libraries I found.

Consider the normal vectorization operation applied to the sparse matrix $G$. Now define a projection matrix $P$ which omits the zero elements to recover what you've denoted as the $g$-vector. So we have
$$eqalign{
g_{vec} &= {rm vec}(G) = P^Tg cr
g &= P,g_{vec} = PP^Tg cr
I &= PP^T cr
}$$
As a concrete example of the projection matrix consider
$$eqalign{
G &= pmatrix{1&4cr 0&0},quad g_{vec}=pmatrix{1cr 0cr 4cr 0},quad g=pmatrix{1cr 4} cr
P &= pmatrix{1&0&0&0cr0&0&1&0},quad P^T=pmatrix{1&0cr 0&0cr 0&1cr 0&0} cr
}$$Define the matrix
$$F = (YG^TLGB - X) = -f$$
Write the energy in terms of $F$ and find its differential in terms of the differential $dg$.
$$eqalign{
{mathcal E} &= |F|^2_F = {rm Tr}(F^TF) = F:F {quadrm Big(Frobenius,productBig)} cr
d{mathcal E}
&= 2F:dF cr
&= 2F:big(Y,dG^T,LGB+YG^TL,dG,Bbig) cr
&= big(2LGBF^TY + 2L^TGY^TFB^Tbig):dG cr
&= 2{,rm vec}big(2LGBF^TY + 2L^TGY^TFB^Tbig):{rm vec}(dG) cr
&= 2Big((Y^TFB^Totimes L),g_{vec} + (BF^TYotimes L^T),g_{vec}Big):dg_{vec} cr
&= 2Big(Y^TFB^Totimes L + BF^TYotimes L^TBig)P^Tg:P^T,dg cr
&= 2PBig(Y^TFB^Totimes L + BF^TYotimes L^TBig)P^Tg:dg cr
}$$
So the gradient is
$$eqalign{
frac{partial {mathcal E}}{partial g} &= 2PBig(Y^TFB^Totimes L+BF^TYotimes L^TBig)P^Tg cr
}$$
Finding the Hessian is not worth the additional effort. Just use a gradient-based method,

e.g. Polak-Ribiere, Barzilai-Borwein, etc.