code: public static DMatrixRMaj[] span(int dimen, int numVectors , Random rand ) {
          if( dimen < numVectors )
              throw new IllegalArgumentException("The number of vectors must be less than or equal to the dimension");
  
          DMatrixRMaj u[] = new DMatrixRMaj[numVectors];
  
          u[0] = RandomMatrices_DDRM.rectangle(dimen,1,-1,1,rand);
          NormOps_DDRM.normalizeF(u[0]);
  
          for( int i = 1; i < numVectors; i++ ) {
  //            System.out.println(" i = "+i);
              DMatrixRMaj a = new DMatrixRMaj(dimen,1);
              DMatrixRMaj r=null;
  
              for( int j = 0; j < i; j++ ) {
  //                System.out.println("j = "+j);
                  if( j == 0 )
                      r = RandomMatrices_DDRM.rectangle(dimen,1,-1,1,rand);
  
                  // find a vector that is normal to vector j
                  // u[i] = (1/2)*(r + Q[j]*r)
                  a.set(r);
                  VectorVectorMult_DDRM.householder(-2.0,u[j],r,a);
                  CommonOps_DDRM.add(r,a,a);
                  CommonOps_DDRM.scale(0.5,a);
  
  //                UtilEjml.print(a);
  
                  DMatrixRMaj t = a;
                  a = r;
                  r = t;
  
                  // normalize it so it doesn't get too small
                  double val = NormOps_DDRM.normF(r);
                  if( val == 0 || Double.isNaN(val) || Double.isInfinite(val))
                      throw new RuntimeException("Failed sanity check");
                  CommonOps_DDRM.divide(r,val);
              }
  
              u[i] = r;
          }
  
          return u;
      }
docs: is there a faster algorithm out there? This one is a bit sluggish