Firstly, solve the initial problem:
$$ \varphi_0 = \underset{\varphi}{\arg\min}\lVert \left( \varphi ^{\top}\varphi -I \right) X \rVert _{F}^{2}, $$
s.t.
$$ \lVert \varphi_0 \rVert _2=1, $$
where $I$ is the $n \times n$ identity matrix, $\varphi_0\in \mathbb{R}^{1\times n}$ is a vector, $X\in \mathbb{R}^{n\times k}$ is a given matrix (dataset).
Secondly, for each $i \in \left\{ 1,2,\cdots, q \right\} $, get the best $\varphi_i\in \mathbb{R}^{1\times n}$ by solving the following problem iteratively:
$$ \varphi_i = \underset{\varphi}{\arg\min}\lVert \left( A_i ^{\top}A_i+\varphi ^{\top}\varphi -I \right) X \rVert _{F}^{2}, $$
s.t.
$$ \lVert \varphi_i \rVert _2=1, $$ where
$$ A_i=\left( \begin{array}{c} \varphi _0\\ \cdots\\ \varphi _{i-1}\\ \end{array} \right) _{i\times n}. $$
The final $A_{q}$ composed of $q$ unit vectors is what I need.
This amounts to an application of the EYM theorem. If $X$ has singular value decomposition $X = U\Sigma V^T$ and $U$ has columns $u_1,\dots,u_n$, then we will find that $$ A_q = \pmatrix{u_1^T\\ \vdots \\ u_q^T}. $$