docs: private Pair<DBIDVar, Double> singleIteration(Relation<V> relationx, Relation<? extends NumberVector> relationy) {
      final int dim = RelationUtil.dimensionality(relationx);
      final int dimy = RelationUtil.dimensionality(relationy);
      assert (dim == 2);
      KNNQuery<V> knnQuery = QueryUtil.getKNNQuery(relationx, getDistanceFunction(), k + 1);
  
      // We need stable indexed DBIDs
      ArrayModifiableDBIDs ids = DBIDUtil.newArray(relationx.getDBIDs());
      // Sort, so we can do a binary search below.
      ids.sort();
  
      // init F,X,Z
      double[][] X = new double[ids.size()][6];
      double[][] F = new double[ids.size()][ids.size()];
      double[][] Y = new double[ids.size()][dimy];
  
      {
        int i = 0;
        for(DBIDIter id = ids.iter(); id.valid(); id.advance(), i++) {
          // Fill the data matrix
          {
            V vec = relationx.get(id);
            double la = vec.doubleValue(0);
            double lo = vec.doubleValue(1);
            X[i][0] = 1.0;
            X[i][1] = la;
            X[i][2] = lo;
            X[i][3] = la * lo;
            X[i][4] = la * la;
            X[i][5] = lo * lo;
          }
  
          {
            final NumberVector vecy = relationy.get(id);
            for(int d = 0; d < dimy; d++) {
              double idy = vecy.doubleValue(d);
              Y[i][d] = idy;
            }
          }
  
          // Fill the neighborhood matrix F:
          {
            KNNList neighbors = knnQuery.getKNNForDBID(id, k + 1);
            ModifiableDBIDs neighborhood = DBIDUtil.newArray(neighbors.size());
            for(DBIDIter neighbor = neighbors.iter(); neighbor.valid(); neighbor.advance()) {
              if(DBIDUtil.equal(id, neighbor)) {
                continue;
              }
              neighborhood.add(neighbor);
            }
            // Weight object itself positively.
            F[i][i] = 1.0;
            final int nweight = -1 / neighborhood.size();
            // We need to find the index positions of the neighbors,
            // unfortunately.
            for(DBIDIter iter = neighborhood.iter(); iter.valid(); iter.advance()) {
              int pos = ids.binarySearch(iter);
              assert (pos >= 0);
              F[pos][i] = nweight;
            }
          }
        }
      }
      // Estimate the parameter beta
      // Common term that we can save recomputing.
      double[][] common = times(transposeTimesTranspose(X, F), F);
      double[][] b = times(inverse(times(common, X)), times(common, Y));
      // Estimate sigma_0 and sigma:
      // sigma_sum_square = sigma_0*sigma_0 + sigma*sigma
      double[][] sigmaMat = times(F, minusEquals(times(X, b), times(F, Y)));
      final double sigma_sum_square = normF(sigmaMat) / (relationx.size() - 6 - 1);
      final double norm = 1 / FastMath.sqrt(sigma_sum_square);
  
      // calculate the absolute values of standard residuals
      double[][] E = timesEquals(times(F, minus(Y, times(X, b))), norm);
  
      DBIDVar worstid = DBIDUtil.newVar();
      double worstscore = Double.NEGATIVE_INFINITY;
      int i = 0;
      for(DBIDIter id = ids.iter(); id.valid(); id.advance(), i++) {
        double err = squareSum(getRow(E, i));
        // double err = Math.abs(E.get(i, 0));
        if(err > worstscore) {
          worstscore = err;
          worstid.set(id);
        }
      }
  
      return new Pair<>(worstid, FastMath.sqrt(worstscore));
    }