In trying to write an integral relation in a discrete manner, I got to an equation of the form
$$MAM x=b$$
where $A$ is a given symmetric matrix, $b$ and $x$ are given vectors and unknown $M$ is a diagonal matrix. The values of all the matrices and vectors can be complex. How can I solve this for $M$?
My attempts
I tried to find patterns that would allow for inverting the order in which the matrices are applied but without a solution. Another thing I tried has been to actually write down explicitly the relation in order to get the system of equations. For the first equation, I got something like
$$m_{11} \left( a_{11}m_{11}x_1 + a_{12}m_{22}x_2 + \cdots + a_{1j}m_{jj}x_j + \cdots + a_{1n}m_{nn} x_n \right) = b_1$$
where $m_{ii}$ are the diagonal entries of $M$, $a_{ij}$ are entries of $A$ and so on for vectors $x$ and $b$, and so on for the other equations. This system of equations is not linear and I have no idea how to solve a non-linear system of equations.
I assume that another way in which the problem can be formulated, although it is not really the same is how can one solve a system of equations like
$$Ax = b/x$$
where $b/x$ should be interpreted as a vector defined by the element-wise division of each pair of points $(b/x)_j = b_j/x_j$.
Best Answer
$ \def\o{{\tt1}}\def\p{\partial} \def\L{\left}\def\R{\right} \def\LR#1{\L(#1\R)} \def\BR#1{\Big(#1\Big)} \def\BBR#1{\Bigg(#1\Bigg)} \def\vec#1{\operatorname{vec}\LR{#1}} \def\diag#1{\operatorname{diag}\LR{#1}} \def\Diag#1{\operatorname{Diag}\LR{#1}} \def\trace#1{\operatorname{Tr}\LR{#1}} \def\qiq{\quad\implies\quad} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\dk#1{\LR{#1_k-#1_{k-\o}}} \def\c#1{\color{red}{#1}} $Let's rename the variable $M\to Y,\,$ then the independent variable for this problem is the vector $y$ $$\eqalign{ Y &= \Diag{y} \quad\iff\quad y=\diag{Y} \\ }$$ Define the vector-valued function $$\eqalign{ f(y) &= \BR{YAYx-b} &= \BR{A\odot{yy^T}}x - b \\ }$$ where $(\odot)$ denotes the elementwise/Hadamard product.
The current problem can be restated as a system of (mildly) nonlinear equations $${f(y) = 0}$$ The simplest route to a solution is probably the Barzilai-Borwein method, since (in this situation) it does not require any gradient calculations.
The basic Barzilia-Borwein method is extremely straightforward
Initialize $$\eqalign{ y_0 &= random \qquad\qquad\qquad\qquad\qquad\quad \\ }$$ First step $$\eqalign{ f_0 &= f(y_0) \\ y_1 &= y_0 - \LR{\frac{0.0\o\,\|y_0\|}{\|f_0\|}}f_0 \qquad\qquad\quad\quad \\ k &= \o \\ }$$ Subsequent steps $$\eqalign{ f_k &= f(y_k) \\ y_{k+\o} &= y_k - \LR{\frac{\dk{y}^T\dk{f}}{\dk{f}^T\dk{f}}}f_k \\ k &= k+\o \\\\ }$$ Stop when the function is nearly zero $\;\|f_k\|\le 10^{-12} $
or when the steplength gets too small $\;\BBR{\frac{\|y_k-y_{k-\o}\|}{\o+\|y_k\|}}\le 10^{-12}$
The introduction of a per iteration difference operator $$\eqalign{ dy_k &= y_k - y_{k-\o} \\ }$$ allows for a more concise description of the algorithm $$\eqalign{ y_{k+\o} &= y_{k} - \LR{\frac{dy_{k}^Tdf_{k}}{df_{k}^Tdf_{k}}}f_{k}, \qquad f_{k+\o} &= f(y_{k+\o}) \\ }$$ Suppressing the $k$-index makes it look less cluttered $$\boxed{\eqalign{ y_+ &= y - \LR{\frac{dy^Tdf}{df^Tdf}}f, \qquad \qquad f_+ &= f(y_+) \\ }}$$