Let $A\in\mathbb R^{m\times n}$, $B\in\mathbb R^{m\times k}$ and suppose $A$ has an SVD. Assuming $\mathcal R(B) \subseteq \mathcal R(A),$ characterize all solutions of the matrix linear equation $$AX = B$$ in terms of the SVD of $A$.
Attempt: We know that SVD of $A$ and $A^+$ are $A=U\Sigma V^T$ and $A^+ = V\Sigma^+U^T$, respectively. Thus, all solutions of the matrix equation $$AX=B$$ are of the form $$X = (V\Sigma^+U^T)B +(I- (V\Sigma^+U^T)(U\Sigma V^T))Y.$$ To verify, simply calculate $AX = B$ with our new $X$. Observe, $$ \begin{equation}\begin{split} AX = (U\Sigma V^T)X &= (U\Sigma V^T)\bigg((V\Sigma^+U^T)B +\big(I- (V\Sigma^+U^T)(U\Sigma V^T)\big)Y\bigg) \\ &= (U\Sigma V^T)V\Sigma^+U^TB +(U\Sigma V^T)(I- V\Sigma^+\Sigma V^T) Y \\ &= U\Sigma \Sigma^+U^TB +(U\Sigma V^T)(\underbrace{I- VIV^T}_{0}) Y \\ &= U\Sigma \Sigma^+U^TB + 0\\ &= B. \end{split}\end{equation} $$
Is this what the question (which seemed trivial) was asking for? Feedback is welcomed!
Note that the characterization of all solutions of $AX=B$ not using the SVD is of the form $$X = A^+B +(I- A^+A)Y$$ for $Y\in\mathbb R^{n\times k}$.
NOTE: '$^+$' denotes Moore-Penrose pseudo-inverse.