Let us assume that matrix $S\in \mathbb{R}^{n\times n}$ is symmetric. We generate a random basis $B = \{b_1,\dots, b_n\}$.
We can assume that entries of each vector $b_i$ are randomly generated by iid normal distribution $𝑁(0,\sigma^2)$. Then we apply the Gram Schmidt process.
Can we say that matrix $S$ is PSD with probability one if for all $i\in [n]$ we have $b_i^T S b_i \ge 0$?
The short answer is no ...
I realized that this is not true by the following example:
Suppose symmetric matrix $S \in \mathbb{R}^{2\times2}$ is full rank. Its SVD decomposition is $S = \lambda_+ v_+ v_+^T + \lambda_- v_- v_-^T$ such that $\lambda_-<0<\lambda_+$ are its eigenvalues and vectors $v$'s are the corresponding eigenvectors.
Suppose $b\in \mathbb{R}^2$ is a unit one random vector with angle $\theta_+$ and $\theta_-$ with the corresponding eigenvectors. Then we have
\begin{equation} \begin{aligned} b^T S b &= \lambda_+ \cos^2 \theta_+ + \lambda_- \cos^2 \theta_- \\ % &= \lambda_+ \cos^2\left(\pi/2 - \theta_-\right) + \lambda_- \cos^2 \theta_-\\ % &=\lambda_+ \sin^2 \theta_- + \lambda_- \cos^2 \theta_-\\ % &= \lambda_+ + (\lambda_- - \lambda_+) \cos^2 \theta_- < 0. \end{aligned} \end{equation}
We can conclude that $\theta_- < \arccos\left(\sqrt{\frac{\lambda_+}{\lambda_+ - \lambda_-}}\right) \le \arccos\left(\sqrt{\frac{\kappa}{\kappa + 1}}\right)$ with $\kappa$ being the condition number of $S$. So, depending on the condition number, the random vector $b$ needs to be tilted towards $v_-$. Thus a uniform distribution does not work.
In a higher dimension, we can divide eigenvalues into two sets of positive and negative ones and use the fact that any linear combination of their corresponding eigenvectors is orthogonal.
Now, the question is how to generate vectors $\pmb{b}$'s having an upper bound on the ratio of max and min absolute value of non-zero eigenvalues?