Problem solving inverse heat-kernel matrix equation

Question

I have the following equation: $H = A + {\frac{1}{2}}L^2 - {\frac{1}{6}}L^3$, with $L = I - A$, where $L$, $A$, and $H$ are all the same sized square matrices, $L$ is symmetric and $I$ denotes the identity matrix. I rewrote the equation to the following: $H = {\frac{1}{6}}A^3 + {\frac{1}{2}}A + {\frac{1}{3}}I$. I want to isolate A but can't seem to find a suitable way. The answer is probably obvious. However, it's been a while since my last linear algebra class. A numerical approach would also be more than welcome if there is no closed-form solution. For anybody interested, the first equation is the 3rd-degree Taylor series approximation of the Heat-Kernel solution for a graph, L denotes the symmetric Laplace matrix, and A denotes the normalized adjacency matrix. I really appreciate any help you can provide.

Edit: I used Matlab's symbolic toolbox to evaluate a symbolic solution. However, this takes up to 20min for a single 9x9 (9 node Graph), which is not feasible, which led me to a faster numeric method. It's not a closed-form solution, but with the fixed-point iteration, one can approximate a solution quite well. I formulated the fixed-point problem as follows: $x_{n+1} = 2H-{\frac{2}{3}}I-{\frac{1}{3}}x_n^3$ with $x$ denoting the normalized adjacency matrix ($A$). I used the following python code to approximate the solution:

import numpy as np

def fixPointIteration(func, x0, iter, tol):
    i = 0
    e = 1

    while(i <= iter and e > tol):
        x = np.asarray(func(x0), dtype=float)

        e = mse(func(x), x)
        i += 1

        x0 = x

    print(f"ran for {i} iterations with error {e}")
    return np.round(x, 10)

def normAdj2Adj(nAdj):
    nAdj[nAdj != 0] = 1
    return nAdj.transpose()

dim = heat.shape[0]

func = lambda A: 2*heat - (2/3)*np.eye(dim) - (1/3)*np.linalg.matrix_power(A, 3)

normalizedAdj = fixPointIteration(func=func, x0=np.eye(dim), iter=1000000, tol=1e-30)

recoveredAdj = normAdj2Adj(normalizedAdj)

I hope this can help somebody.

Answer 1

The solution matrix $A$ can be found in a straight forward manner.

The cubic matrix equation is

$$ H = \frac{1}{6} A^3 + \frac{1}{2} A + \dfrac{1}{3} I $$

Since $A$ (the unknown matrix) is symmetric, it can be diagonalized, and written as $ A = R D R^T $, where $D$ is a diagonal matrix and $R$ is an orthogonal (rotation) matrix.

Therefore the right hand side, becomes

$ R \left( \frac{1}{6} D^3 + \frac{1}{2} D + \frac{1}{3} I \right) R^T $

It follows from this that matrix $H$ (which is symmetric) must have the same matrix $R$ as the right hand side.

Thus our problem reduces to solving $n$ scalar cubic equations. If $D_H$ is the diagonal matrix associated with $H$, then

$ D_H = \frac{1}{6} D^3 + \frac{1}{2} D + \frac{1}{3} I $

Taking the diagonal elements of both sides one by one, a cubic equation in $D_{ii}$ can be constructed, and solved. Since this is a cubic equation, we have at least one real root, for $i = 1, \dots ,n $

Once all the diagonal elements of matrix $D$ have been computed, matrix $A$ is given by

$ A = R D R^T $

I've tested the above method and verified that it gives the correct solutions for the matrix $A$.