Edit: I made a mistake in my numerical experiments, as numpy.linalg.svd returns the "canonical" full SVD, where $ \mathbf{U}, \mathbf{V} \in \mathbb{R}^{d \times d}$ are both unitary, $ \mathbf{S} \in \mathbb{R}^{d \times d} $, although many of the approximated singular values are very small, i.e. on the order of $ 10^{-18} $, and so can be discarded. Using only the first $ r \approx \operatorname{rank}\mathbf{C} $ columns of $ \mathbf{V} $, $ r $ computed by numpy.linalg.matrix_rank using default parameters, numpy.allclose returns False, a more sensible result. In general, numerical experiments indicate that $ \mathbf{VV}^\top \ne \mathbf{I} $ when $ \operatorname{rank}\mathbf{V} < d $.
Let $ \mathbf{C} \in \mathbb{R}^{d \times d} $ be such that $ \mathbf{C} \succeq \mathbf{0} $, i.e. $ \mathbf{C} $ positive semidefinite, where $ \operatorname{rank}\mathbf{C} = r < d $. Suppose we compute the singular value decomposition [1] of $ \mathbf{C} $, i.e. we have $ \mathbf{C} = \mathbf{USV}^\top $, where $ \mathbf{U}, \mathbf{V} \in \mathbb{R}^{d \times r} $ with orthonormal columns, $ \mathbf{S} \in \mathbb{R}^{r \times r} $ diagonal containing only the nonzero singular values of $ \mathbf{C} $. How "close" is $ \mathbf{VV}^\top $ to $ \mathbf{I} $? For now I have tentatively concluded that given a linear system $ \mathbf{Cx} = \mathbf{b} $, denoting the Moore-Penrose pseudoinverse of $ \mathbf{C} $ as $ \mathbf{C}^+ $, that it is ok to write $ \mathbf{w} \approx \mathbf{C}^+\mathbf{b} $ when $ \mathbf{C} $ is singular. Of course, if $ \mathbf{C} \succ \mathbf{0} $, i.e. $ \mathbf{C} $ positive definite, $ \mathbf{w} = \mathbf{C}^{-1}\mathbf{b} $.
Any thoughts and help would be greatly appreciated. Thanks!
[1] I use the definition as stated in A.5.4 of Boyd and Vandenberghe's Convex Optimization.