The Schmidt decomposition is just a confusing red herring here.
Given a basis $\lvert \psi_i\rangle$ of our Hilbert space, we may expand any vector in that basis, and hence a normalized state is $\lvert \psi\rangle = \sum c_i\lvert \psi_i\rangle$ with $\sum_i \lvert c_i\rvert^2 = 1$. So its associated density matrix is
$$ \rho_\psi = \lvert \psi\rangle\langle \psi\rvert = \sum_{i,j} c_ic^\ast_j \lvert \psi_i\rangle\langle \psi_i\rvert = \sum_{i,j} \rho_{ij}\lvert \psi_i\rangle\langle \psi_j\rvert. \tag{1}$$
You seem to be confused because this "looks" as if there is more freedom here in the $\rho_{ij}$ than in the $\lambda_i$ when we write a generic mixed state as
$$ \rho = \sum_i \lambda_i \lvert \phi_i\rangle \langle \phi_i\rvert,\tag{2}$$
in terms of its eigenstates $\lvert \phi_i\rangle$ - after all there are $n^2$ of the $\rho_{ij}$ and only $n$ of the $\lambda_i$.
There are a few different ways to note that eq. (2) defines a more general class of operators than eq. (1):
Eq. (1) is an expansion of an arbitrary $\rho_\psi$ in a fixed basis $\lvert \psi_i\rangle$, something that works for all pure state matrices, while eq. (2) is specifically in the eigenbasis of $\rho$ - some other $\rho'$ would not have the same form in that same basis in general.
There isn't actually more independent parameters in eq. (1) because there are only as much of the $c_i$ as there are of the $\lambda_i$ in eq. (2), and while the $\lambda_i$ are constrained by $\sum_i \lambda_i = 1$, the $c_i$ are constrained by $\sum_i \lvert c_i\rvert^2 = 1$, so in both cases you essentially have $n-1$ independent parameters, but in eq. (2) you can additionally choose the basis differently and this gives you different $\rho$s (see the first point again), so there is more freedom in eq. (2) than in eq. (1).
A very "manifest" difference is simply that the pure state matrix fulfills the purity condition $\rho_\psi^2 = \rho_\psi$ by construction, while the mixed state matrix does not - it does not follow from eq. (2) and $\sum_i \lambda_i = 1$, while it does follow from eq. (1) and $\sum_i \lvert c_i\rvert^2 = 1$, so that alone shows you the second equation is more general.
Your claim already fails in the simplest non-trivial case, namely for three qubits.
This follows from a (the?) classic result in the theory of multipartite entanglement -- Three qubits can be entangled in two inequivalent ways by Dür, Vidal, and Cirac.
They show that for three qubits, there are two classes of genuine tripartite entangled states, the W class and the GHZ class. Given states $\lvert \phi_\mathrm{W}\rangle$ and $\lvert \phi_\mathrm{GHZ}\rangle$, from the two classes, it is impossible to convert between them using SLOCC, i.e., it is impossible to write
$$
\lvert\phi_\mathrm{W}\rangle\propto (A\otimes B\otimes C)\lvert\phi_\mathrm{GHZ}\rangle
$$
(or vice versa) for any $A$, $B$, and $C$. In particular, this rules out the possibility of unitaries doing the job, as you are asking.
What remains is to show that there are states $\lvert\phi_\mathrm{W}\rangle$ and $\lvert \phi_\mathrm{GHZ}\rangle$ with identical Schmidt coefficients in every bipartition. This is equivalent to demanding that their single-qubit reduced states have the same spectra.
Using the results of the paper of Dür, Vidal, and Cirac, it is now easy to find such states (and indeed, they really should exist -- it would be rather disappointing if the two classes were distinguished by the spectra of their reduced density matrices). For instance, you can choose
$$
\lvert \phi_\mathrm{W}\rangle =
\sqrt{\gamma}\lvert 001\rangle +
\sqrt{\gamma}\lvert 010\rangle +
\sqrt{\gamma}\lvert 100\rangle +
\sqrt{1-3\gamma}\lvert 000\rangle
$$
(cf. Eq. (20) in the paper), and
$$
\lvert \phi_\mathrm{GHZ}\rangle
\propto
\lvert 0\rangle \lvert 0\rangle \lvert 0\rangle +
\lvert \theta\rangle\lvert \theta\rangle\lvert \theta\rangle
$$
(cf. Eq. (15)),
with $\lvert\theta\rangle=\sqrt{\mu}\lvert 0\rangle + \sqrt{1-\mu}\lvert 1\rangle$.
You can now easily check that for the spectra $(\lambda,1-\lambda)$ ($\lambda\le1/2$) of the single-qubit reduced states (which are all equal by symmetry)
- for $\lvert\phi_\mathrm{W}\rangle$, all values $0<\lambda\le1/3$ can be obtained by varying $0<\gamma\le1/3$
- for $\lvert\phi_\mathrm{GHZ}\rangle$, all values $0<\lambda\le1/2$ can be obtained by varying $0\le\mu< 1$.
Thus, you can easily find values $\mu$ and $\gamma$ where the reduced states have identical spectra, and thus the Schmidt spectra in all bipartitions are equal, yet the states cannnot be converted into each other by local unitaries (or even SLOCC).
Best Answer
Concretely, given an object/tensor with any number of indices, you can choose any bipartition of the indices and perform the corresponding SVD. More precisely, this means to think of the object as an operator acting between two linear spaces, and perform its SVD as you would normally.
For example, your $A$ is an operator acting in a bipartite space, and therefore has $4$ indices. Write these as $A=\sum_{ijk\ell}A_{ijk\ell} |ij\rangle\!\langle k\ell|$. The "standard SVD" corresponds to thinking of this object as the linear operator $$|k\ell\rangle\mapsto \sum_{ij} A_{ijk\ell} |ij\rangle,$$ and doing the SVD of the corresponding matrix. However, you could also consider the linear operator $$|j\ell\rangle\mapsto \sum_{ik} A_{ijk\ell} |ik\rangle.$$ If you perform the SVD of the corresponding matrix (which is obtained from the other one transposing some layers), you get a decomposition of the form $$A_{ijk\ell} = \sum_m s_m u^m_{ik} \,\bar v^m_{j\ell},$$ for some $s_m\ge0$ and collections of orthonormal vectors $\{u^m\}_m,\{v^m\}_m$. If you use this decomposition in the original formula for $A$, you get $$A = \sum_{ijk\ell} \sum_m s_m u^m_{ik}\,\bar v^m_{j\ell} |ij\rangle\!\langle k\ell| = \sum_m s_m B_m\otimes C_m,$$ where $$B_m\equiv \sum_{ik} u^m_{ik} |i\rangle\!\langle k|, \qquad C_m\equiv \sum_{ik} \bar v^m_{j\ell} |j\rangle\!\langle \ell|.$$
Note that the orthonormality of the vectors $u^m,v^m$, given by the SVD, translate into the orthonormality of the operators $B_m,C_m$ with respect to the Hilbert-Schmidt inner product: $$\operatorname{Tr}(B_m^\dagger B_n)=\operatorname{Tr}(C^\dagger_m C_n)=\delta_{mn}.$$
Now, if you want to do this operator decomposition in practice, given an operator expressed as a matrix, you need to "reshape" the matrix to obtain a four-index tensor object, then convert this into a new matrix, and compute the SVD of that. If you have access to Mathematica, here's a quick snippet to take a matrix and return the corresponding operator Schmidt decomposition:
and here's an example of what you get in a few simple cases: