docs: public static CovarianceMatrix make(Relation<? extends NumberVector> relation) {
      int dim = RelationUtil.dimensionality(relation);
      CovarianceMatrix c = new CovarianceMatrix(dim);
      double[] mean = c.mean;
      int count = 0;
      // Compute mean first:
      for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
        NumberVector vec = relation.get(iditer);
        for(int i = 0; i < dim; i++) {
          mean[i] += vec.doubleValue(i);
        }
        count++;
      }
      if(count == 0) {
        return c;
      }
      // Normalize mean
      for(int i = 0; i < dim; i++) {
        mean[i] /= count;
      }
      // Compute covariances second
      // Two-pass approach is numerically okay and fast, when possible.
      double[] tmp = c.nmea; // Scratch space
      double[][] elems = c.elements;
      for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
        NumberVector vec = relation.get(iditer);
        for(int i = 0; i < dim; i++) {
          tmp[i] = vec.doubleValue(i) - mean[i];
        }
        for(int i = 0; i < dim; i++) {
          for(int j = i; j < dim; j++) {
            elems[i][j] += tmp[i] * tmp[j];
          }
        }
      }
      // Restore symmetry.
      for(int i = 0; i < dim; i++) {
        for(int j = i + 1; j < dim; j++) {
          elems[j][i] = elems[i][j];
        }
      }
      c.wsum = count;
      return c;
    }