stat/mds: reduce work in TorgersonScaling and fix k return

The fix for k papers over inaccuracies in near zero Eigenvalues. So this
may want further work there.

stat/mds/mds.go

@@ -62,18 +62,25 @@ func TorgersonScaling(dst *mat.Dense, eigdst []float64, dis mat.Symmetric) (k in
 	vals := ed.Values(nil)
 	reverse(vals, dst.RawMatrix())
 	copy(eigdst, vals)
+
+	// FIXME(kortschak): The Eigenvalues near zero are probably zero.
+	// Without this, the final column of the eurodist test case is 0
+	// and the returned k differs from the result from R's cmdscale.
+	const tol = 1e-10
+
 	for i, v := range vals {
-		if v < 0 {
+		if v < tol {
 			vals[i] = 0
 			continue
 		}
 		k = i + 1
 		vals[i] = math.Sqrt(v)
 	}
 
-	// TODO(kortschak): Make use of the knowledge
-	// of k to avoid doing unnecessary work.
-	dst.Mul(dst, mat.NewDiagDense(len(vals), vals))
+	var tmp mat.Dense
+	tmp.Mul(dst, mat.NewDiagonalRect(n, k, vals[:k]))
+	dst = dst.Slice(0, n, 0, k).(*mat.Dense)
+	dst.Copy(&tmp)
 
 	return k, dst, eigdst
 }

stat/mds/mds_test.go

@@ -71,7 +71,7 @@ var torgersonScalingTests = []struct {
 			3927, 2868, 1616, 1786, 2196, 1403, 650, 2068, 3886, 949, 1500, 3231, 2108, 3188, 2428, 2187, 1754, 1827, 2707, 0, 2105,
 			1991, 1802, 1175, 1381, 1588, 937, 1455, 1019, 2974, 1155, 1205, 2937, 1157, 2409, 1363, 898, 428, 1249, 1209, 2105, 0,
 		}),
-		wantK: 12,
+		wantK: 11,
 		want: mat.NewDense(21, 11, []float64{
 			2290.274680, 1798.80293, 53.79314, -103.826958, -156.955115, 54.755434, -47.6768205, 1.241284, -14.893196, -6.366664, 4.818373,
 			-825.382790, 546.81148, -113.85842, 84.585831, 291.440759, -33.046236, -74.5267190, 3.766233, 225.620420, -21.270973, 22.050484,
@@ -114,14 +114,12 @@ func TestTorgersonScaling(t *testing.T) {
 			t.Error("unexpected scaling failure")
 			continue
 		}
-		_, wc := test.want.Dims()
-		n := test.dis.Symmetric()
 		if gotK != test.wantK {
 			t.Errorf("unexpected k for test %d: got:%d want:%d", i, gotK, test.wantK)
 		}
-		if !mat.EqualApprox(colAbs{got.Slice(0, n, 0, wc)}, colAbs{test.want}, 1e-6) {
+		if !mat.EqualApprox(colAbs{got}, colAbs{test.want}, 1e-6) {
 			t.Errorf("unexpected result for test %d:\ngot:\n%.4f\nwant:\n%.4f",
-				i, mat.Formatted(got.Slice(0, n, 0, wc)), mat.Formatted(test.want))
+				i, mat.Formatted(got), mat.Formatted(test.want))
 		}
 		if !floats.EqualApprox(gotVals, test.wantVals, 1e-10) {
 			t.Errorf("unexpected Eigenvalues for test %d:\ngot: %.12e\nwant:%.12e", i, gotVals, test.wantVals)