public double update(double[] x) {
        if (x.length != n) {
            throw new IllegalArgumentException(String.format("Invalid input vector size: %d, expected: %d", x.length, n));
        }

        projection.mv(x, y);

        for (int j = 0; j < p; j++) {
            for (int i = 0; i < n; i++) {
                double delta = x[i];
                for (int l = 0; l <= j; l++) {
                    delta -= projection.get(l, i) * y[l];
                }
                projection.add(j, i, r.apply(t) * y[j] * delta);

                if (Double.isInfinite(projection.get(j, i))) {
                    throw new IllegalStateException("GHA lost convergence. Lower learning rate?");
                }
            }
        }

        t++;
        projection.mv(x, y);
        projection.tv(y, wy);
        return MathEx.squaredDistance(x, wy);
    }
@Test
    public void test() {
        System.out.println("GHA");

        int k = 3;
        double[] mu = MathEx.colMeans(USArrests.x);
        Matrix cov = Matrix.of(MathEx.cov(USArrests.x));
        for (int i = 0; i < USArrests.x.length; i++) {
           MathEx.sub(USArrests.x[i], mu);
        }

        TimeFunction r = TimeFunction.constant(0.00001);
        GHA gha = new GHA(4, k, r);
        for (int iter = 1, t = 0; iter <= 1000; iter++) {
            double error = 0.0;
            for (int i = 0; i < USArrests.x.length; i++, t++) {
                error += gha.update(USArrests.x[i]);
            }
            error /= USArrests.x.length;

            if (iter % 100 == 0) {
                System.out.format("Iter %3d, Error = %.5g%n", iter, error);
            }
        }

        Matrix p = gha.projection;
        Matrix t = p.ata();
        System.out.println(t);

        Matrix s = p.mm(cov).mt(p);
        double[] ev = new double[k];
        System.out.println("Relative error of eigenvalues:");
        for (int i = 0; i < k; i++) {
            ev[i] = Math.abs(eigenvalues[i] - s.get(i, i)) / eigenvalues[i];
            System.out.format("%.4f ", ev[i]);
        }
        System.out.println();

        for (int i = 0; i < k; i++) {
            assertTrue(ev[i] < 0.1);
        }

        double[][] lt = MathEx.transpose(loadings);
        double[] evdot = new double[k];
        double[][] pa = p.toArray();
        System.out.println("Dot products of learned eigenvectors to true eigenvectors:");
        for (int i = 0; i < k; i++) {
            evdot[i] = MathEx.dot(lt[i], pa[i]);
            System.out.format("%.4f ", evdot[i]);
        }
        System.out.println();

        for (int i = 0; i < k; i++) {
            assertTrue(Math.abs(1.0- Math.abs(evdot[i])) < 0.1);
        }
    }