I'm not an expert on these kind of "optimization" problems, but here is what I can say.
[EDIT.] I think Part 2. of this "algorithm" is in fact quite useless, because, if you already have a basis for your subspace $V$ you'd better apply Part 3. which I'm adding now.
Part 1. First of all, if your subspace $V \subset \mathbb{R}^N$ was given by a set of cartesian equations,
$$
AX = 0 \ ,
$$
then, getting the reduced row-echelon form of matrix $A$,
$$
\begin{pmatrix}
1 & 0 & \dots & 0 & a^1_{r+1} & \dots & a^1_n \\
0 & 1 & \dots & 0 & a^2_{r+1} & \dots & a^2_n \\
\dots \\
0 & 0 & \dots & 1 & a^r_{r+1} & \dots & a^r_n \\
0 & 0 & \dots & 0 & 0 & \dots & 0 \\
\dots \\
0 & 0 & \dots & 0 & 0 & \dots & 0
\end{pmatrix}
$$
would allow you to obtain some of the unknowns as function of the others. Here $r$ is the rank of $A$. Namely,
\begin{eqnarray*}
x_1 &=& - (a^1_{r+1} x_{r+1} + \dots + a^1_{n}x_n) \\
x_2 &=& - (a^2_{r+1} x_{r+1} + \dots + a^2_{n}x_n ) \\
\dots \\
x_r &=& - (a^r_{r+1} x_{r+1} + \dots + a^r_{n}x_n )
\end{eqnarray*}
So, you would obtain a basis for your vector subspace with lots of zeros, because any vector of it could be written as
$$
(x_1, \dots , x_r, x_{r+1}, \dots , x_n) = x_{r+1}(-a^1_{r+1}, \dots, -a^r_{r+1}, 1, 0, \dots , 0) + x_2(-a^2_1, \dots, -a^2_n, 0, 1, 0, \dots , 0) +\dots + x_{n} (-a^r_1, \dots , -a^r_n, 0, \dots , 0, 1) \ .
$$
That is,
\begin{eqnarray*}
u_1 &=& (-a^1_{r+1}, \dots, -a^r_{r+1}, 1, 0, \dots , 0) \\
u_2 &=& (-a^2_1, \dots, -a^2_n, 0, 1, 0, \dots , 0) \\
\dots && \\
u_d &=& (-a^r_1, \dots , -a^r_n, 0, \dots , 0, 1) \ ,
\end{eqnarray*}
with $d = n-r$, would be a basis for your vector subspace with a pretty amount of zeros everywhere.
Part 2. But you don't start with a system of cartesian equations, do you? -Instead, you already have a bais for your $V$. Well, then you should obtain first a system of linear equations for it. For instance, like this: if $v_1, \dots , v_d$ is your basis, then any vector $u\in V$ can be written as
$$
u = \lambda_1 v_1 + \dots + \lambda_d v_d \ .
$$
(Hence $d$ is the dimension of $V$.) If
$$
v_1 =
\begin{pmatrix}
a^1_1 \\
\vdots \\
a^n_1
\end{pmatrix}
\ , \dots ,
v_d =
\begin{pmatrix}
a^1_d \\
\vdots \\
a^n_d
\end{pmatrix}
$$
were the coordinates of your basis vectors, then this would mean that the $\lambda_i$ up there of your
$$
u =
\begin{pmatrix}
x_1 \\
\vdots \\
x_n
\end{pmatrix}
$$
would be the solution of this linear system of equations:
$$
\begin{pmatrix}
a^1_1 & \dots & a^1_d \\
\vdots & & \vdots \\
a^n_1 & \dots & a^n_d
\end{pmatrix}
\begin{pmatrix}
\lambda_1 \\
\vdots \\
\lambda_d
\end{pmatrix}
=
\begin{pmatrix}
x_1 \\
\vdots \\
x_n
\end{pmatrix} \ .
$$
Ok, so now reduce to row-echelon form this system,
$$
\begin{pmatrix}
* & \dots & * & | & g_1(x_1, \dots , x_n) \\
\vdots & &\vdots & | & \\
* & \dots & * & | & g_d(x_1, \dots , x_n) \\
0 & \dots & 0 & | & g_{d+1}(x_1, \dots , x_n) \\
\vdots & &\vdots & | & \\
0 & \dots & 0 & | & g_n(x_1, \dots , x_n) \\
\end{pmatrix}
$$
and you'll get a linear system of equations for your subspace. Namely,
$$
g_{d+1}(x_1, \dots , x_n)= 0 \ , \dots , \ g_n(x_1, \dots , x_n) = 0 \ ,
$$
to which you can apply part 1.
Example. Let's say you have $V \subset \mathbb{R}^3$ given by a basis like
$$
v_1 = (1,-1,0) , v_2 = (1,1, 0) \ .
$$
So, clearly $V$ is the $xy$-plane and we could safely produce a "simple" basis with no computations at all, just taking
\begin{eqnarray*}
u_1 &=& (1,0,0) \\
u_2 &=& (0,1,0)
\end{eqnarray*}
But let's see how our "algorithm" works, nevertheless. First we use the second part to produce a system of linear equations for $V$:
$$
\begin{pmatrix}
1 & 1 & | & x \\
-1 & 1 & | & y \\
0 & 0 & | & z
\end{pmatrix}
$$
Adding the first row to the second, you get
$$
\begin{pmatrix}
1 & 1 & | & x \\
0 & 2 & | & x +y \\
0 & 0 & | & z
\end{pmatrix}
$$
Which is already in row-echelon form, so functions $g_i$ are in this case
$$
g_1(x,y,z) = x \ , \ g_2(x,y,z) = x + y \quad \text{and} \quad g_3(x,y,z) = z \ .
$$
The only one that matters for us is the last one: $z = 0$. Now, we apply part 1 and get from here that all vectors of $V$ have the following form
$$
(x,y,0) = x (1,0,0) + y (0,1,0) \ ,
$$
for arbitrary values of $x,y$. So the "simple" basis of $V$ is, indeed,
$$
u_1 = (1,0,0) \quad \text{and} \quad u_2 = (0,1,0) \ .
$$
Part 3. If you have a basis $v_1, \dots , v_d$ for $V$ with coordinates as in Part 2., you just put together all your vectors as the columns of a matrix
\begin{pmatrix}
a^1_1 & \dots & a^1_d \\
\vdots & & \vdots \\
a^n_1 & \dots & a^n_d
\end{pmatrix}
and then reduce it to its reduced row-echelon form, but by means of elementary column transformations, instead of row ones. Since elementary transformations are the same as linear combinations and they preserve the rank of the matrix, you won't change the subspace generated by $v_1, \dots , v_d$, that is, $V$. And in the end, you'll get a basis with plenty of zeros too. Something like
$$
\begin{pmatrix}
1 & 0 & \dots & 0 \\
0 & 1 & \dots & 0 \\
\vdots & \vdots & \ddots & 0 \\
0 & 0 & \dots & 1 \\
b^d_1 & b^d_2 & \dots & b^d_d \\
\vdots & \vdots & \ddots & \vdots \\
b^n_1 & b^n_2 & \dots & b^n_d
\end{pmatrix}
$$
Applying this new "algorithm" to the same example, this time we get:
$$
\begin{pmatrix}
1 & 1 \\
-1 & 1 \\
0 & 0
\end{pmatrix}
\longrightarrow
\begin{pmatrix}
1 & 0 \\
-1 & 2 \\
0 & 0
\end{pmatrix}
\longrightarrow
\begin{pmatrix}
1 & 0 \\
-1 & 1 \\
0 & 0
\end{pmatrix}
\longrightarrow
\begin{pmatrix}
1 & 0 \\
0 & 1 \\
0 & 0
\end{pmatrix}
$$
Here the elementary column transformations have been the following:
- Substract the first column from the second one.
- Divide the second column by $2$.
- Add the second column to the first one.
