Assuming that we have a matrix $A\in \mathbb{R}^{m\times n}$ stored in a quantum superposition, i.e. $$|A\rangle= \frac{1}{\|A\|_F}\sum_{i,j=0}^{n-1}{a_{ij}}|i,j\rangle$$ and a vector $b\in \mathbb{R}^{n\times 1}$ stored in $$|b\rangle= \frac{1}{\|b\|}\sum_{i=0}^{n-1}{b_{i}}|i\rangle.$$
How can I find a unitary, such that
$$U|b\rangle\approx|Ab\rangle?$$
I do not know what the state-of-the-art for this problem is, but here is my attempt at it.
First, I doubt whether the required unitary $U$ can be found for every matrix $A$. The most glaring issue occurs when $b\in\ker A$. In this case, $|Ab\rangle$ is ill-defined, so it is not clear what $U|b\rangle$ should be approximately equal to. A slightly more subtle issue concerns the preservation of the inner product. Consider for example $A=\begin{pmatrix}1 & 0 \\ 0 & 2\end{pmatrix}$ and $b_1=(1, 0)^T$ and $b_2=(1, 1)^T$. Then $|b_1\rangle = |0\rangle$, $|b_2\rangle = |+\rangle$ and $|Ab_1\rangle = |0\rangle$, $|Ab_2\rangle = (|0\rangle + 2|1\rangle)/\sqrt{5}$. However, for any unitary $U$ we have
$$ \langle b_1|U^\dagger U|b_2\rangle = \langle b_1|b_2\rangle = \frac{1}{\sqrt{2}} $$
which is hard to regard as approximately equal to
$$ \langle Ab_1|Ab_2\rangle = \frac{1}{\sqrt{5}}. $$
Such issues may be resolved by using the context in which the problem arises to place additional restrictions on $A$ and $b$. Alternatively, we can widen the search for the required transformation beyond the unitary operations. The latter is what we do below.
We prepare the state $|Ab\rangle$ by teleporting the state $|b\rangle$ through the state $|A\rangle$ and post-selecting on a measurement outcome.
Suppose we have three subsystems $K$, $L$ and $M$ each with Hilbert space of dimension $n$. First, prepare the subsystem $K$ in the state $|b\rangle=\frac{1}{\|b\|}\sum_{k=0}^{n-1}|k\rangle_K$ and the subsystems $L$ and $M$ in the joint state $|A\rangle=\frac{1}{\|A\|_F}\sum_{m,l=0}^{n-1}a_{ml}|l\rangle_L|m\rangle_M$. Now, measure subsystems $K$ and $L$ in the Fourier basis
$$ |\Psi_{u,v}\rangle_{KL} = \frac{1}{\sqrt{n}} (X^u Z^v \otimes I) \sum_{j=0}^{n-1} |j\rangle_K|j\rangle_L\tag1 $$
where $X$ and $Z$ are the generalization of the Pauli operators to $n$-dimensional Hilbert space known as the shift and clock matrices and given by $X|j\rangle = |j + 1 \mod n\rangle$ and $Z|j\rangle = e^{\frac{2\pi ij}{n}}|j\rangle$. See this answer for a good description of qudit teleportation.
Each possible measurement outcome is described by a pair of numbers $u, v \in \{0, 1, \dots, n-1\}$. Suppose that the measurement outcome is $u=0, v=0$. Then the subsystem $M$ is in the state
$$ \begin{align} _{KL}\langle\Psi_{0,0}|b\rangle_K|A\rangle_{LM} &=\left(\frac{1}{\sqrt{n}}\sum_{k=0}^{n-1}{}_K\langle k| {}_L\langle k|\right) |b\rangle_K|A\rangle_{LM}\\ &=\frac{1}{\sqrt{n}}\sum_{k=0}^{n-1}{}_K\langle k|b\rangle_K\, {}_L\langle k|A\rangle_{LM} \\ &=\frac{1}{\sqrt{n}}\sum_{k=0}^{n-1}\frac{b_k}{\|b\|} \sum_{l=0}^{n-1} \frac{a_{lk}}{\|A\|_F}|l\rangle_M \\ &=\frac{1}{\sqrt{n}\|b\|\|A\|_F}\sum_{l=0}^{n-1} \sum_{k=0}^{n-1}a_{lk}b_k|l\rangle_M \\ &=\frac{1}{\sqrt{n}\|b\|\|A\|_F}\sum_{l=0}^{n-1} (Ab)_l|l\rangle_M \\ &=\frac{\|Ab\|}{\sqrt{n}\|b\|\|A\|_F}|Ab\rangle_M \end{align}\tag2 $$
which after normalization becomes $|Ab\rangle$ as promised. For other measurement outcomes the subsystem $M$ is not in the state $|Ab\rangle$ due to the operators $X$ and $Z$ in $(1)$. In some applications it may be possible to boost success probability by applying a correction based on the measurement outcome $u, v$ to obtain $|Ab\rangle$. However, this part of the protocol depends on properties of $A$ and $b$.
It is interesting to consider what happens in certain special cases. For example, consider the situation where $b \in \ker A$. In this case, the state $|Ab\rangle$ is ill-defined and from $(2)$ we see that the probability of obtaining the required $u=0,v=0$ measurement outcome is zero.
Another interesting case occurs when $A=aa^T$ is a rank one matrix and $a^Tb\ne 0$. Then $|A\rangle_{LM}=|a\rangle_L|a\rangle_M$ is a product state and nothing can be teleported through it. However in this case
$$ |aa^Tb\rangle = \frac{1}{\|aa^Tb\|}\sum_{k=0}^{n-1} a_k a^Tb |k\rangle = \frac{a^Tb}{|a^Tb|}|a\rangle \equiv |a\rangle $$
and the initial state is the appropriate output state up to the unobservable global phase.