Minimizing error for $A$ in $Ax=b$ if $x$ and $b$ are matrices There was an answer give for how to solve $Ax=b$ if $x$ and $b$ are known here:
Solving for $A$ in $Ax = b$
But I was wondering is there a way to solve it if there is no exact solution and you want to minimize error, similar to how you solve $Ax=b$ for $x$ if using least squares if there is no exact solution.
Additionally, say $X$ and $B$ are matrices that correspond to a bunch of $Ax=b$ equations. Are there known methods for minimizing error in that case?
Thanks for any help.
 A: This is an interesting question in my opinion. It is a baby step towards the low-rank Hessian update formulas that are used in quasi-Newton methods (SR1, DFP, etc).
Here is one way to go about it:
First, let us frame it as the following optimization problem:
$$\boxed{\min_A \frac{1}{2}\|AX - B\|_\text{Fro}^2 + \frac{\alpha}{2}\|A\|_\text{Fro}^2}$$
where $\|\cdot\|_\text{Fro}$ is the Frobenius matrix norm (square root of the sum of squares of all matrix entries), and $\alpha$ is a regularization parameter.
It is straightforward to see that $\|M\|_\text{Fro}^2 = \operatorname{trace}\left[M^TM\right]$, where $\operatorname{trace}\left[\cdot\right]$ is the matrix trace (sum of all diagonal entries). Thus the objective function may be written as
$$\frac{1}{2}\operatorname{trace}\left[\left(AX - B\right)^T\left(AX - B\right)\right] + \frac{\alpha}{2}\operatorname{trace}\left[A^TA\right].$$
The solution to the optimization problem is the matrix $A$ such that the derivative of the objective function is zero if $A$ is perturbed in any direction $H$. Using the linearity and cyclic properties of the trace, we have
\begin{align*}
0 &= \frac{d}{dA}\left(\frac{1}{2}\operatorname{trace}\left[\left(AX - B\right)^T\left(AX - B\right)\right] + \frac{\alpha}{2}\operatorname{trace}\left[A^TA\right]\right) \cdot H \\
&= \operatorname{trace}\left[\left(HX\right)^T\left(AX - B\right)\right]+\alpha \operatorname{trace}\left[H^TA\right]  \\
&= \operatorname{trace}\left[H^T\left(\left(AX - B\right)X^T + \alpha A\right)\right]
\end{align*}
Since this must hold for all perturbations $H$, we must have
$$0 = \left(AX - B\right)X^T + \alpha A$$
or
$$A\left(XX^T + \alpha I\right) = BX^T$$
which is a linear system that you may solve to find the desired matrix $A$.
What's more, using the Woodbury formula we have
$$\left(XX^T + \alpha I\right)^{-1} = \frac{1}{\alpha}\left( I - X\left(\alpha I + X^TX\right)^{-1}X^T\right)$$
and therefore
$$\boxed{A = \frac{1}{\alpha}BX^T - \frac{1}{\alpha}BX^T X\left(\alpha I + X^TX\right)^{-1}X^T .}$$
Here is some python/numpy code that tests the formula:
import numpy as np

N = 10
r = 3

X = np.random.randn(N,r)
B = np.random.randn(N,r)

a = 1e-6
print('a=', a)
A = (1./a) * B @ X.T @ (np.eye(N) - X @ np.linalg.inv(a*np.eye(r) + X.T @ X) @ X.T)

rel_res = np.linalg.norm(A @ X - B) / np.linalg.norm(B)
print('rel_res=', rel_res)

norm_A = np.linalg.norm(A)
print('norm_A=', norm_A)

and here is the output of running that code:
a= 1e-06
rel_res= 2.3549592377658484e-07
norm_A= 2.24920744100694

