Given the matrices ${\bf A}, {\bf B},{\bf C}$, how to solve $\bf A\Lambda B = C$ for diagonal $\bf \Lambda$?

Question

Given the fat matrix ${\bf A} \in \mathbb{C}^{M \times Q}$ (where $M < Q$), the tall matrix ${\bf B} \in \mathbb{C}^{Q \times N}$ (where $Q > N$) and ${\bf C} \in \mathbb{C}^{M \times N}$, how to solve the following linear matrix equation for diagonal matrix ${\bf \Lambda} \in \mathbb{C}^{Q \times Q}$?

$$ {\bf A} {\bf \Lambda} {\bf B} = {\bf C} $$

Answer 1

$ \def\LR#1{\left(#1\right)} \def\op#1{\operatorname{#1}} \def\diag#1{\op{diag}\LR{#1}} \def\Diag#1{\op{Diag}\LR{#1}} \def\trace#1{\op{Tr}\LR{#1}} \def\frob#1{\left\| #1 \right\|_F} \def\qiq{\quad\implies\quad} \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\c#1{\color{red}{#1}} \def\CLR#1{\c{\LR{#1}}} \def\fracLR#1#2{\LR{\frac{#1}{#2}}} \def\gradLR#1#2{\LR{\grad{#1}{#2}}} $First, introduce the extremely useful Frobenius $(:)$ product $$\eqalign{ F:G &= \sum_{i=1}^m\sum_{j=1}^n F_{ij}G_{ij} \;=\; \trace{F^TG} \\ F^*:F &= \frob{F}^2 \qquad \{ {\rm Frobenius\;norm} \}\\ F:G &= G:F \;=\; G^T:F^T \\ H:\LR{FG} &= \LR{HG^T}:F \;=\; \LR{F^TH}:G \\ }$$ and a diagonalization identity involving the Hadamard $(\odot)$ product $$ \diag{A\cdot\Diag{x}\cdot B} \;\equiv\; \LR{B^T\odot A}x $$ Then approach this as a least squares problem $$\eqalign{ &M = A\cdot\Diag{x}\cdot B \;-\; C \\ &\min_x\: \frob{M}^2 }$$ Calculate the gradient of the least squares objective function wrt the $x$ vector $$\eqalign{ \phi &= {M^*:M} \\ d\phi &= M^*:\c{dM} \\ &= M^*:\CLR{A\cdot\Diag{dx}\cdot B} \\ &= \LR{A^TM^*B^T}:\Diag{dx} \\ &= \diag{A^TM^*B^T}:dx \\ &= \diag{A^TA^*\cdot\Diag{x^*}\cdot B^*B^T - A^TC^*B^T}:dx \\ \grad{\phi}{x} &= \LR{BB^H\odot A^TA^*}x^* - \diag{A^TC^*B^T} \\ }$$ Conjugate the gradient, set it to zero and solve for the optimal $x$ $$\eqalign{ x &= \LR{B^*B^T\odot A^HA}^{\boldsymbol +}\diag{A^HCB^H} \\ {\large\Lambda} &= \Diag{x} \\ }$$ where
$\quad M^+\:$ denotes the pseudoinverse of $M$
$\quad M^T,M^*,M^H\:$ denote the transpose, complex and hermitian conjugates
$\quad\Diag{}\:$ creates a diagonal matrix from its vector argument
$\quad\diag{}\;$ returns the main diagonal of its matrix argument as a column vector

Answer 2

By comparing the $(m,n)$ element of matrices $AB$ and $C$, note that

$$C_{m,n} = \sum_{i=1}^Q A_{m,i} \Lambda_{i,i} B_{i,n}$$

This gives $M\cdot N$ linear equations in $Q$ unknowns ($\Lambda_{1,1}, \Lambda_{2,2}, \ldots, \Lambda_{Q,Q}$), if $\Lambda$ is really a diagonal matrix.

Answer 3

First, notice that this is essentially just a set of linear equations, since $$A\operatorname{diag}(x+y)B = A\operatorname{diag}(x)B + A\operatorname{diag}(y)B$$

Secondly, notice that there are $M×N$ many equations, but $Q$ variables. So depending on whether $M×N ≷ Q$, the system will be over- or underdetermined. When $M×N > Q$, the system generally will not be solveable and the best you can hope for is a minimal residual solution like least-squares. Letting $X= \operatorname{diag}(x) = (δ_{pq}x_q)_{pq}$ we have:

$$\begin{aligned} AXB &= \Big(∑_{pq} A_{mp}X_{pq}B_{qn}\Big)_{mn} \\&= \Big(∑_{q} A_{mq}B_{qn}x_q\Big)_{mn} \\&= \Big(A_{mq}B_{qn}\Big)_{mn, q} ⋅ x \end{aligned}$$

After which we can flatten the first two axes and use any off-the shelf solver. Example in numpy:

import numpy as np

M, N, Q = 2, 3, 6
A = np.random.randn(M, Q)
B = np.random.randn(Q, N)
C = np.random.randn(M, N)


def lstsq(A, B, C):
    """Returns least-squares solution for A⋅diag(x)⋅B = C."""
    M = np.einsum("mq, qn -> mnq", A, B)
    M = M.reshape(-1, Q)  # flatten
    C = C.reshape(-1)  # flatten
    return np.linalg.lstsq(M, C, rcond=None)


x, *_ = lstsq(A, B, C)
diff = np.einsum("mq, q, qn -> mn", A, x, B) - C
print(diff)

Answer 4

Assume that there is a unique algebraic solution of the equation $$ C = A\Lambda B $$ Multiplying by $\{X,Y\}$ does not affect the equality, but creates square-shaped matrices $$ XCY = X(A\Lambda B)Y $$ The main diagonals of each side are also equal $$\eqalign{ {\rm diag}(XCY) &= {\rm diag}(XA\Lambda BY) &= \Big((BY)^T \odot(XA)\Big)\,{\rm diag}(\Lambda) \\ }$$ This simple linear system is easy to solve $$ {\rm diag}(\Lambda) = \Big((BY)^T\odot(XA)\Big)^{-1}\,{\rm diag}(XCY) $$ Any two matrices $\{X,Y\}$ of the proper dimensions can be used, but the choice will affect the invertibility of the linear system. Some obvious candidates are $$\eqalign{ X \in \{A^H,A^T\} \qquad Y \in \{B^H,B^T\} \\ }$$ If the original equation does not have a solution, then you must settle for an approximate least-squares solution. In that case, Greg's analysis indicates that $\,\{X=A^H,\:Y=B^H\}\,$ are the appropriate choices.