docs: public int correlationDistance(PCAFilteredResult pca1, PCAFilteredResult pca2, int dimensionality) {
      // TODO: Can we delay copying the matrixes?
      // pca of rv1
      double[][] v1t = copy(pca1.getEigenvectors());
      double[][] v1t_strong = pca1.getStrongEigenvectors();
      int lambda1 = pca1.getCorrelationDimension();
  
      // pca of rv2
      double[][] v2t = copy(pca2.getEigenvectors());
      double[][] v2t_strong = pca2.getStrongEigenvectors();
      int lambda2 = pca2.getCorrelationDimension();
  
      // for all strong eigenvectors of rv2
      double[][] m1_czech = pca1.dissimilarityMatrix();
      for(int i = 0; i < v2t_strong.length; i++) {
        double[] v2_i = v2t_strong[i];
        // check, if distance of v2_i to the space of rv1 > delta
        // (i.e., if v2_i spans up a new dimension)
        double distsq = squareSum(v2_i) - transposeTimesTimes(v2_i, m1_czech, v2_i);
  
        // if so, insert v2_i into v1 and adjust v1
        // and compute m1_czech new, increase lambda1
        if(lambda1 < dimensionality && distsq > deltasq) {
          adjust(v1t, v2_i, lambda1++);
          // TODO: make this incremental?
          double[] e1_czech_d = new double[v1t.length];
          Arrays.fill(e1_czech_d, 0, lambda1, 1);
          m1_czech = transposeDiagonalTimes(v1t, e1_czech_d, v1t);
        }
      }
  
      // for all strong eigenvectors of rv1
      double[][] m2_czech = pca2.dissimilarityMatrix();
      for(int i = 0; i < v1t_strong.length; i++) {
        double[] v1_i = v1t_strong[i];
        // check, if distance of v1_i to the space of rv2 > delta
        // (i.e., if v1_i spans up a new dimension)
        double distsq = squareSum(v1_i) - transposeTimesTimes(v1_i, m2_czech, v1_i);
  
        // if so, insert v1_i into v2 and adjust v2
        // and compute m2_czech new , increase lambda2
        if(lambda2 < dimensionality && distsq > deltasq) {
          adjust(v2t, v1_i, lambda2++);
          // TODO: make this incremental?
          double[] e2_czech_d = new double[v1t.length];
          Arrays.fill(e2_czech_d, 0, lambda2, 1);
          m2_czech = transposeDiagonalTimes(v2t, e2_czech_d, v2t);
        }
      }
      return Math.max(lambda1, lambda2);
    }