Quote:
Originally Posted by wgardner
Yeah, that. For a real positive semi-definite symmetric matrix (like A'A for any A), the SVD is something like U D U' where U is orthogonal and D is diagonal. In our case, A'A is not full rank, so the last diagonal value of D is 0.
|
With rounding error, it will not be exactly zero. For the 2015 MISJO event, it is 3.43E-16. So your code needs to use a threshold.
Quote:
|
Using the method in the link above, the pseudo-inverse is computed as U E U' where E is diagonal with elements Ei = 1/Di except where Di=0
|
use abs(Di)<threshold instead of Di==0
Quote:
|
One other formulation for CCWMOA
|
We need a better acronym. That's too awkward.
Quote:
... would just be:
if we have T teams, have T-1 unknown values C1, C2, ..., C(T-1) and set CT = -Sum(C1, C2,... C(T-1)) in all of the equations (thus enforcing that all T values of Ci are zero mean). Then we only have T-1 equations with T-1 unknowns and everything is full rank. This is just another way of saying we want to find the values of C1, C2, ... CT that minimize the prediction error subject to the constraint that the resulting set of Ci values have zero mean.
|
I ran the numbers for all 117 events in 2015 and found that the following two computational methods yield virtually identical results for min L2 norm of b-Ax:
Method 1
1a) [U,S,V] = svd(A)
1b) Replace the diagonal elements of S with their reciprocals, except when abs(Sjj)<threshold (I used 1e-4 for threshold), in which case make Sjj zero.
1c) compute x = V*S*(U'*b)
Method 2
2a) N = A'*A
2b) d= A'*b
2c) compute x = N\d
..... (Octave mldivide notation)
2d) compute m = mean(x)
2e) subtract m from each element of x
Notice Method 1 factors A, not A'A, resulting in less rounding error.