I need to solve the following equation for $X$ with $d$-by-$d$ matrices $A,B,C,D$ and Hadamard product $\odot$

$$AXB + X\odot C = D$$

Vectorizing all terms gives a solution with $O(d^6)$ complexity, which is intractable since $d\approx 1000$ in my application. Is there something I can do to get an estimate in $O(d^3)$ time?

To add more information about the structure, A,B,C are moment matrices. Specifically, for random variables $X,Y$ they are

$$A_{ij}=E[X_iX_j]$$ $$B_{ij}=E[Y_iY_j]$$ $$C_{ij}=E[X_iY_j]$$ $$E_t[D_{ij}]=C_{ij}$$

To give even more background, this comes up in problem of speeding up neural network training which for a single layer can be viewed as the following problem

$$\text{minimize}_{W} E[(x'Wy)^2]$$

Here $x$ and $y$ are random variables with shape $(d,1)$, the gradient is $E[x_iy_j]$ and the curvature is $E[x_i x_j y_k y_l]$. The goal is to use curvature information to obtain Newton-like correction to a noisy estimate of the gradient. Applying Isserlis theorem we can approximate curvature rank-4 tensor in terms of rank-2 covariance tensors which leads to matrices $A,B,C$ above. Substituting small sample estimate of gradient into $D$ and large sample estimate of curvature into $A,B,C$, then solving for $X$, gives us a preconditioned gradient step.

Incidentally, using the same sample to estimate curvature and gradient, gives slightly simplified problem:

$$AXB + X\odot C = C$$

$A$ and $B$ are known to be ill-conditioned -- $n$-th eigenvalue is approximately $1/n$ and majority of eigenvalues are numerically 0, possibly 90% of all eigenvalues. Coincidentally, this ill-conditioning is what allows neural networks to generalize. This implies that $C$ is also singular.

  • 1
    $\begingroup$ Can you clarify what $\odot$ means? $\endgroup$ Sep 28 '19 at 3:26
  • 2
    $\begingroup$ I think what to do depends on the data $A, B, C$. If $A,B$ are invertible and if $\|A^{-1}\| \|B^{-1}\|\|C\|<1$ (say for the Frobenius norms) then $X\mapsto X + A^{-1}(X\odot C)B^{-1}$ is invertible by a Neumann series, with the usual estimate for the remainder. On the other hand, if all coefficients of $C$ are non-zero, and $\|A \| \|B \|\|C^{-\odot}\|<1$ (here $C^{-\odot} $ denotes the Hadamard inverse) one can invert instead $X\mapsto X + (AXB)\odot C^{-\odot} $. In both cases the above conditions on the norms are of course stronger than needed for the convergence. $\endgroup$ Sep 28 '19 at 7:45
  • 2
    $\begingroup$ So I think a relevant point is: are your $A,B,C,D$ totally generic $d\times d$ matrices, or do they have any special feature that may address to a convenient method? $\endgroup$ Sep 28 '19 at 7:56
  • 1
    $\begingroup$ I would try some iterative method like conjugated gradient possible with a preconditioner but the choice depends on $A,B,C$ $\endgroup$
    – user35593
    Sep 28 '19 at 10:58
  • 1
    $\begingroup$ Cross-posted on Computational Science. $\endgroup$ Sep 30 '19 at 11:25

I am going to reduce your problem to another form. In a truncated version, this modified problem has been discussed here and it seems there is no apparent efficient solution à la Bartels–Stewart. However, it still might be valuable to have an alternative approach at hand, especially because such equations are being actively discussed (see T. Damm, Direct methods and ADI-preconditioned Krylov subspace methods for generalized Lyapunov equations, Numer. Linear Algebra Appl. 15, no. 9 (2008): 853-871) and because iterative approaches based on the standard Lyapunov solvers do exist.

Let us apply the vec-operation to the both sides of your equation using the commutativity of Hadamard product

$$ \mathrm{vec}(AXB)+\mathrm{vec}(C\odot X)=\mathrm{vec}(D). $$

Now we use the following properties $$ \text{vec}(ABC)=(C^\mathrm{T}\otimes A)\text{vec}(B) $$ and $$ \text{vec}(A\odot B)=\text{vec}(A)\odot\text{vec}(B). $$ Thus $$ (B^\mathrm{T}\otimes A)\text{vec}(X)+\mathrm{diag}\!\left[\mathrm{vec} (C)\right]\mathrm{vec} (X)=\mathrm{vec}(D). $$ Let us introduce a new matrix $U$: $$ U=(B^\mathrm{T}\otimes A)+\mathrm{diag}\!\left[\mathrm{vec} (C)\right]. $$ Our goal is now to write it in the form of a Kronecker product. $$ U=\sum_i \sigma_i V_i^\mathrm{T}\otimes W_i. $$ This is known as the nearest Kronecker product problem. Using SVD decomposition on a permuted version of $U$ as Van Loan (J. Comp. Appl. Math, 123 (2000) 85) proposed, the factors $V_i$ and $W_i$ and the singular values $\sigma_i$ can be determined, and the original equation can be written as

$$ \sum_i \sigma_i\left(V_i^\mathrm{T}\otimes W_i\right)\mathrm{vec} (X)=\mathrm{vec}(D), $$ which is equivalent to $$ \sum_i \sigma_i W_i X V_i=D. $$ Truncating the sum to just 2 terms, a standard Lyapunov equation is obtained.

  • 1
    $\begingroup$ This would just be a preconditioner though, right? $\endgroup$
    – Nick Alger
    Oct 1 '19 at 3:34
  • $\begingroup$ @NickAlger Why do you say so? $\endgroup$
    – yarchik
    Oct 1 '19 at 7:48
  • 2
    $\begingroup$ It seems to me that only $U \approx V^T\otimes W$, and they cannot be made equal in general (unless I am misunderstanding something). For example imagine $A=B=I$, and $C$ is a matrix filled with iid gaussian random entries. Then $U$ is diagonal, so $V$ and $W$ must be diagonal. But then equality is not possible, by a degree of freedom counting argument ($2d$ vs $d^2$). $\endgroup$
    – Nick Alger
    Oct 1 '19 at 8:14
  • $\begingroup$ @NickAlger Yes, you are right. $\endgroup$
    – yarchik
    Oct 1 '19 at 8:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.