I was wondering if there is exists an (efficient) way to implement the state
$$|\Psi\rangle=\left(\sum_i^N c_i|{\psi_i}\rangle\right)=\left(\sum_i^N c_i\hat U_i\right)|{0.....0}\rangle$$
where $\hat{U}_i$ are known unitaries and the states $|{\psi_i}\rangle, |{\psi_j}\rangle$ are most likely not orthogonal, that is $\langle\psi_j|{\psi_i}\rangle\neq \delta_{ij}$.
You might be able to use something like the linear combination of unitaries protocol from (Childs, 2017). Summarizing Lemmas 6-7 in that paper, you will need any unitary $V$ that can prepare the following state on $m = \lceil \log_2 N \rceil$ qubits:
$$ V|0^m\rangle = \frac{1}{\sqrt{c}} \sum_{i=1}^N \sqrt{c_i} |i\rangle \tag{1} $$ where $c = \sum_i c_i$, and you need the ability to implement the unitary \begin{equation} U = \sum_{i=1}^N |i\rangle \langle i| \otimes U_i \tag{2} \end{equation} which you can think of as a controlled gate generalized to qudits - see answers to this question for more discussion on such gates. Then assuming that your unitaries $U_i$ each act on $n$ qubits you can show that
$$ (V^\dagger \otimes I) U (V \otimes I)|0^m\rangle |0^n\rangle = \frac{1}{c} |0^m\rangle |\Psi\rangle + |\Phi^\perp\rangle \tag{3} $$
where $|\Psi\rangle$ is the state you're trying to prepare. The $|\Phi^\perp\rangle$ appears when you implement $V^\dagger$ - since we only specify the first column in defining $V$, applying $V^\dagger$ gives us a "garbage" output associated with every other basis state in the first register besides $|0^m\rangle$.
Once you've prepared the state in Eq. $(3)$, you can just measure the first $m$ qubits and with probability $(\lVert |\Psi\rangle\rVert/c)^2 $ you will observe $0^m$, which means you've prepared $|\Psi\rangle$ in the second register. Alternatively, you can gain a square-root speedup by first performing amplitude amplification on the first term of the state given in Eq. $(3)$.
