The climate in k numbers
A city’s climate in the article on principal component analysis is described by 24 numbers: the mean temperature and total precipitation of each month. Less needs to be stored. It is allowed to replace each city by a point on some k-dimensional plane in 24-dimensional space, not necessarily through zero and not necessarily along the axes, and to restore the city by that point. What is the smallest possible mean squared miss?
Formally: given n points x1, …, xn in a space of dimension d and a number k. For a k-dimensional plane P the miss of a point is the squared distance from it to P. Find the smallest value over all k-dimensional planes of the mean miss
$$E_k = \min_P \frac{1}{n}\sum_{i=1}^{n} \operatorname{dist}^2(x_i, P).$$A 0-dimensional plane is a single point. Tests three to six are the climate of 469 cities from the article: temperature in tenths of a degree, precipitation in millimetres.
Input
The first line holds n, d, k (1 ≤ n ≤ 2000, 1 ≤ d ≤ 30, 0 ≤ k ≤ d). Then n lines of d integers each, of absolute value at most 10⁵.
Output
One number: Ek, with an absolute or relative error of 10⁻⁶.
Examples
Example 1
4 2 1 0 1 2 3 4 4 6 7
0.091180794
Example 2
3 3 0 1 2 3 4 0 6 7 8 2
20.444444444
Limits
- Time · Python
- 6000 ms
- Time · Go
- 1000 ms
- Memory
- 256 MB
- Judging
- 10 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
Start with k = 0: which point best replaces the cloud, and what is the mean miss then? How is it expressed in terms of the covariance matrix?
2Hint 2
If the k columns with the greatest variance are kept, the miss is the sum of the variances of the rest. But the plane may run at an angle to the axes. In which axes is the covariance matrix diagonal?
3The idea
The best plane passes through the centre of mass of the cloud: shifting towards it reduces the miss of any plane, as for k = 0, where the best point is the mean and the miss is the sum of the variances, the trace of the covariance matrix S = (1/n)·Σ(xi − x̄)(xi − x̄)T. The sum of squared distances to the points of a plane splits by Pythagoras into a “shadow” on the plane and a “miss” across it. The more shadow, the smaller the miss, and the sum of shadows on a plane with an orthonormal basis v1, …, vk equals ΣvjTSvj. The largest value of this sum is the sum of the k largest eigenvalues of S, attained on the principal components. So
$$E_k = \lambda_{k+1} + \lambda_{k+2} + \dots + \lambda_d,$$the sum of the discarded eigenvalues. This is the statement of Eckart and Young (1936) in the form for clouds of points.
4Walkthrough
Centre the columns, build the d × d matrix S and find all its eigenvalues. For a symmetric 30 × 30 matrix the method of Jacobi rotations is convenient: each rotation in the plane of two coordinates zeroes a pair of off-diagonal elements, the sum of their squares falls quickly from pass to pass, and the eigenvalues remain on the diagonal. A dozen passes are enough. The power method with deflation also works, but when eigenvalues are close it converges slowly, and here all of them are needed.
The sum of the discarded values is better added directly than as the trace minus the sum of the kept ones: in the second case a small answer comes from subtracting two large numbers. On the cities’ climate two components leave a mean miss several times smaller than the two best column axes; this is seen by comparing with the wrong solution below.
4+Solution
5Complexity
O(n·d²) for the covariance matrix and O(d³) per Jacobi pass, with a handful to a few dozen passes.
6Common mistakes
- Keeping columns rather than directions. Choosing the k columns with the greatest variance and discarding the rest is natural, but it is just one of the planes, parallel to the axes. Temperatures of neighbouring months are strongly related, and a plane at an angle wins many times over.
- Not centring. The eigenvalues of the matrix ΣxixiT/n without subtracting the mean describe planes through zero, and the best plane need not pass through zero.
- Losing a small answer. With k = d the answer is exactly 0; in the seventh test the cloud lies in a 4-dimensional plane, and the answer is zero too, computed up to rounding.