I want to find $x\in\mathbb{R^n}$, where $x$ minimizes
$$f(x) = \frac{1}{2}x^TA x + b^Tx$$
There are no constraints. I do know that $A$ is symmetric positive semidefinite, and $f(x) \ge 0$ for all $x$. If $A$ is invertible, then the problem is easy, but I am interested in a general solution that handles the case where $A$ is singular. If the solution is not unique, that's ok. I just need any solution, but preferably a sparse one.
Through some reading, it looks like the pseudoinverse obtained via SVD would give a solution, but I haven't seen much theory about that. Would that work? And if it does, can I get a proof? Or is there another way to do this?
Edit: forgot to mention that $A$ is symmetric
As I'm sure you know, the solutions are all vectors $x$ satisfying the optimality conditions: $$\nabla f(x) = 0 \quad\Longleftrightarrow\quad Ax + b = 0$$ It's important to note that if $b\not\in\mathop{\textrm{Range}}(A)$, the objective is unbounded below, and there is no minimum.
But if $b\in\mathop{\textrm{Range}}(A)$, one solution is $x=-A^{\dagger} b$, where $A^{\dagger}$ represents the Moore-Penrose pseudoinverse. If $A=U_1\Sigma_1 U_1^T$ is the "economy" Schur decomposition of $A$—that is, where $\Sigma_1$ only contains the non-zero eigenvalues values—then $A^\dagger = U_1\Sigma_1^{-1} U_1^T$.
A practical approach to solving this problem might be to compute $x = -A^{\dagger} b$ per the technique above, examine the size of the residual $Ax+b$, and see if it is "small" enough to be numerically considered zero---that is, whether it indicates that $b\in\mathop{\textrm{Range}}(A)$ to within numerical tolerances.
If you want to know all of the solutions, then we need to a bit more work. Specifically, $$x\in\left\{y - A^{\dagger} b \,\middle|\, Ay=0 \right\}$$ so we need a way to characterize the null space of $A$. Let's go back to the Schur decomposition, but this time compute the whole thing: $$A = \begin{bmatrix} U_1 & U_2 \end{bmatrix} \begin{bmatrix} \Sigma_1 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} U_1^T \\ U_2^T \end{bmatrix}$$ It turns out that $U_1$ is an orthonormal basis for that null space. So we have $$x\in\left\{U_2 z_1 - U_1 \Sigma_1^{-1} U_1^T b \,\middle|\, z_1\in\mathbb{R}^p \right\}$$ where $p$ is the number of zero eigenvalues.