nx_matrix_eig_gate.nx source

1// nx_matrix_eig_gate.nx -- AM1: DOES THE SYMMETRIC EIGENSOLVER ACTUALLY SOLVE, AND CAN IT REFUSE? 2// 3// The done-rule was pre-declared on /compare/appliedmath before a line of the solver was written: 4// "reconstructing the matrix from its factors returns the original within the declared bound, the 5// eigenvalues of a known symmetric matrix match their exact values, and a singular matrix is reported 6// singular rather than inverted" 7// so these teeth are the rule being met rather than a rule written to fit what the code happens to do. 8// 9// EVERY TOOTH IS AN INVARIANT OF THE EIGENPROBLEM, NOT A STORED ANSWER. Trace equals the sum of the 10// eigenvalues; a diagonal matrix has its own diagonal as its spectrum; eigenvectors of a symmetric matrix 11// are orthonormal; and the residual max|AV - VL| is the bound the reference libraries publish. A stored 12// golden vector would record that the author agreed with the code on one afternoon; these hold for any 13// correct implementation and fail for any wrong one. 14// 15// license_tier: ORIGINAL expect_exit: 0 16import "syscalls.nx" 17import "nx_matrix.nx" 18import "nx_gate_verdict.nx" 19 20const EG_Q: i64 = 16384 21// DECLARED TOLERANCE, in Q14 units, for values that pass through rotations. Twelve sweeps of fixed-point 22// plane rotations accumulate rounding; 64 units is 0.4 percent of one, and it is stated here rather than 23// discovered by a reader, so nobody can mistake approximate for exact. 24const EG_TOL: i64 = 64 25// The residual bound the solver is required to meet on these fixtures. Published as a NUMBER the gate 26// emits, so a later regression shows up as the figure moving rather than as a boolean flipping. 27const EG_RESID_MAX: i64 = 96 28 29func eg_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v } 30 31func main() -> i64 { 32 let ctr: *i64 = gv_ctr() 33 gv_head("=== NX-MATRIX-EIG-GATE -- symmetric eigenvalues by cyclic Jacobi, in fixed point, with a published residual ===" as *u8) 34 35 // ---- A KNOWN SPECTRUM: [[2,1],[1,2]] has eigenvalues exactly 3 and 1 --------------------------- 36 let a: *Matrix = nx_matrix_alloc(2, 2) 37 nx_matrix_set(a, 0, 0, 2 * EG_Q); nx_matrix_set(a, 0, 1, 1 * EG_Q) 38 nx_matrix_set(a, 1, 0, 1 * EG_Q); nx_matrix_set(a, 1, 1, 2 * EG_Q) 39 let v: *Matrix = nx_matrix_alloc(2, 2) 40 let ev: *i64 = (sys_mmap(16)) as *i64 41 let rc: i64 = nx_matrix_eig_sym(a, ev, v, 12) 42 43 gv_check_eq("fixture-reached-the-condition: the solver ACCEPTS a symmetric matrix and returns OK" as *u8, 44 rc, NXM_EIG_OK, ctr) 45 46 var hi: i64 = ev[0] 47 var lo: i64 = ev[1] 48 if lo > hi { hi = ev[1]; lo = ev[0] } 49 gv_check_near("the-LARGER-eigenvalue-of-a-known-2x2-is-exactly-3-within-the-declared-tolerance" as *u8, 50 hi, 3 * EG_Q, EG_TOL, ctr) 51 gv_check_near("the-SMALLER-eigenvalue-of-a-known-2x2-is-exactly-1-within-the-declared-tolerance" as *u8, 52 lo, 1 * EG_Q, EG_TOL, ctr) 53 54 // TRACE IS INVARIANT UNDER SIMILARITY: the eigenvalues must sum to the trace, whatever they are. 55 // This holds for ANY symmetric matrix, so it cannot be satisfied by a solver that returns the diagonal. 56 gv_check_near("TRACE-INVARIANCE-the-eigenvalues-sum-to-the-trace" as *u8, 57 ev[0] + ev[1], nx_matrix_get(a, 0, 0) + nx_matrix_get(a, 1, 1), EG_TOL, ctr) 58 59 // ---- THE RESIDUAL, WHICH IS THE WHOLE POINT OF PUBLISHING A BOUND ------------------------------ 60 let resid: i64 = nx_matrix_residual(a, v, ev) 61 gv_check("RECONSTRUCTION-A-times-V-equals-V-times-Lambda-within-the-published-residual-bound" as *u8, 62 (resid <= EG_RESID_MAX) as i64, ctr) 63 64 // NEG-CONTROL FOR THE RESIDUAL ITSELF. A residual function that returned zero for everything would 65 // make the tooth above vacuous, so a deliberately corrupted eigenvalue must drive it far up. 66 let bad: *i64 = (sys_mmap(16)) as *i64 67 bad[0] = ev[0] + (EG_Q / 2) 68 bad[1] = ev[1] 69 let resid_bad: i64 = nx_matrix_residual(a, v, bad) 70 gv_check("neg-control-a-CORRUPTED-eigenvalue-drives-the-residual-far-above-the-bound" as *u8, 71 (resid_bad > EG_RESID_MAX * 4) as i64, ctr) 72 73 // ---- EIGENVECTORS OF A SYMMETRIC MATRIX ARE ORTHONORMAL ---------------------------------------- 74 var col0: i64 = nxm_fmul(nx_matrix_get(v, 0, 0), nx_matrix_get(v, 0, 0)) 75 + nxm_fmul(nx_matrix_get(v, 1, 0), nx_matrix_get(v, 1, 0)) 76 gv_check_near("each-EIGENVECTOR-is-a-UNIT-vector-so-the-rotations-preserved-length" as *u8, 77 col0, EG_Q, EG_TOL, ctr) 78 let dot: i64 = nxm_fmul(nx_matrix_get(v, 0, 0), nx_matrix_get(v, 0, 1)) 79 + nxm_fmul(nx_matrix_get(v, 1, 0), nx_matrix_get(v, 1, 1)) 80 gv_check_near("distinct-EIGENVECTORS-are-ORTHOGONAL-which-only-a-real-diagonalisation-gives" as *u8, 81 dot, 0, EG_TOL, ctr) 82 83 // ---- A DIAGONAL MATRIX IS ITS OWN SPECTRUM ------------------------------------------------------ 84 let d: *Matrix = nx_matrix_alloc(3, 3) 85 nx_matrix_fill(d, 0) 86 nx_matrix_set(d, 0, 0, 5 * EG_Q) 87 nx_matrix_set(d, 1, 1, 2 * EG_Q) 88 nx_matrix_set(d, 2, 2, 9 * EG_Q) 89 let vd: *Matrix = nx_matrix_alloc(3, 3) 90 let evd: *i64 = (sys_mmap(24)) as *i64 91 let rcd: i64 = nx_matrix_eig_sym(d, evd, vd, 12) 92 gv_check_eq("a-DIAGONAL-3x3-is-accepted" as *u8, rcd, NXM_EIG_OK, ctr) 93 gv_check_near("a-DIAGONAL-matrix-returns-its-own-diagonal-as-its-spectrum-summed-invariantly" as *u8, 94 evd[0] + evd[1] + evd[2], 16 * EG_Q, EG_TOL, ctr) 95 96 // ---- A SINGULAR MATRIX IS REPORTED SINGULAR, NOT INVERTED --------------------------------------- 97 // [[1,1],[1,1]] has eigenvalues 2 and 0. The ZERO eigenvalue is the report of singularity, and it is 98 // the thing a caller must see before trying to invert. 99 let s: *Matrix = nx_matrix_alloc(2, 2) 100 nx_matrix_set(s, 0, 0, 1 * EG_Q); nx_matrix_set(s, 0, 1, 1 * EG_Q) 101 nx_matrix_set(s, 1, 0, 1 * EG_Q); nx_matrix_set(s, 1, 1, 1 * EG_Q) 102 let vs: *Matrix = nx_matrix_alloc(2, 2) 103 let evs: *i64 = (sys_mmap(16)) as *i64 104 let rcs: i64 = nx_matrix_eig_sym(s, evs, vs, 12) 105 gv_check_eq("a-SINGULAR-matrix-is-accepted-rather-than-rejected" as *u8, rcs, NXM_EIG_OK, ctr) 106 var small: i64 = evs[0] 107 if eg_abs(evs[1]) < eg_abs(evs[0]) { small = evs[1] } 108 gv_check_near("a-SINGULAR-matrix-REPORTS-a-zero-eigenvalue-so-a-caller-learns-not-to-invert-it" as *u8, 109 small, 0, EG_TOL, ctr) 110 gv_check_eq("and-the-determinant-agrees-that-it-is-singular" as *u8, 111 nx_matrix_det_2x2(s), 0, ctr) 112 113 // ---- REFUSALS, EACH NAMING ITS OWN RULE --------------------------------------------------------- 114 let asym: *Matrix = nx_matrix_alloc(2, 2) 115 nx_matrix_set(asym, 0, 0, 1 * EG_Q); nx_matrix_set(asym, 0, 1, 5 * EG_Q) 116 nx_matrix_set(asym, 1, 0, 0); nx_matrix_set(asym, 1, 1, 1 * EG_Q) 117 let va: *Matrix = nx_matrix_alloc(2, 2) 118 let eva: *i64 = (sys_mmap(16)) as *i64 119 gv_check_eq("neg-control-an-ASYMMETRIC-matrix-is-REFUSED-by-name-because-Jacobi-would-return-plausible-nonsense" as *u8, 120 nx_matrix_eig_sym(asym, eva, va, 12), NXM_EIG_ERR_ASYMMETRIC, ctr) 121 122 let rect: *Matrix = nx_matrix_alloc(2, 3) 123 gv_check_eq("neg-control-a-NON-SQUARE-matrix-is-REFUSED-by-its-own-named-rule" as *u8, 124 nx_matrix_eig_sym(rect, eva, va, 12), NXM_EIG_ERR_NOT_SQUARE, ctr) 125 126 gv_check_eq("neg-control-the-symmetry-predicate-SEPARATES-a-symmetric-matrix-from-an-asymmetric-one" as *u8, 127 nx_matrix_is_symmetric(a) - nx_matrix_is_symmetric(asym), 1, ctr) 128 129 // ---- SVD: THE SECOND HALF OF AM1 --------------------------------------------------------------- 130 // The teeth are INVARIANTS of the factorisation, so they hold for any correct SVD and fail for any 131 // wrong one: the squared singular values sum to the Frobenius norm, their product is the absolute 132 // determinant, and U*diag(sv)*V^T reconstructs the original. 133 let sa: *Matrix = nx_matrix_alloc(2, 2) 134 nx_matrix_set(sa, 0, 0, 2 * EG_Q); nx_matrix_set(sa, 0, 1, 1 * EG_Q) 135 nx_matrix_set(sa, 1, 0, 1 * EG_Q); nx_matrix_set(sa, 1, 1, 2 * EG_Q) 136 let su: *Matrix = nx_matrix_alloc(2, 2) 137 let svv: *Matrix = nx_matrix_alloc(2, 2) 138 let svals: *i64 = (sys_mmap(16)) as *i64 139 let rcv: i64 = nx_matrix_svd(sa, su, svals, svv, 12) 140 gv_check_eq("fixture-reached-the-condition: the SVD accepts a general matrix and returns OK" as *u8, 141 rcv, NXM_SVD_OK, ctr) 142 143 gv_check("SINGULAR-VALUES-come-back-in-DESCENDING-order-which-rank-truncation-depends-on" as *u8, 144 (svals[0] >= svals[1]) as i64, ctr) 145 gv_check("SINGULAR-VALUES-are-NON-NEGATIVE-by-definition" as *u8, 146 (svals[1] >= 0) as i64, ctr) 147 148 // FROBENIUS IDENTITY: the squared singular values sum to the sum of the squared entries. True of any 149 // correct SVD of any matrix, and impossible to satisfy by returning the diagonal or the eigenvalues. 150 let frob_sv: i64 = nxm_fmul(svals[0], svals[0]) + nxm_fmul(svals[1], svals[1]) 151 var frob_a: i64 = 0 152 var fr: i64 = 0 153 while fr < 2 { 154 var fc: i64 = 0 155 while fc < 2 { frob_a = frob_a + nxm_fmul(nx_matrix_get(sa, fr, fc), nx_matrix_get(sa, fr, fc)); fc = fc + 1 } 156 fr = fr + 1 157 } 158 gv_check_near("FROBENIUS-IDENTITY-the-squared-singular-values-sum-to-the-sum-of-squared-entries" as *u8, 159 frob_sv, frob_a, EG_TOL * 4, ctr) 160 161 // DETERMINANT IDENTITY: the product of the singular values is the absolute determinant. 162 gv_check_near("DETERMINANT-IDENTITY-the-product-of-the-singular-values-is-the-absolute-determinant" as *u8, 163 nxm_fmul(svals[0], svals[1]), 3 * EG_Q, EG_TOL * 4, ctr) 164 165 let sresid: i64 = nx_matrix_svd_residual(sa, su, svals, svv) 166 gv_check("RECONSTRUCTION-U-times-Sigma-times-V-transpose-returns-the-original-within-the-bound" as *u8, 167 (sresid <= EG_RESID_MAX) as i64, ctr) 168 169 let sbad: *i64 = (sys_mmap(16)) as *i64 170 sbad[0] = svals[0] + (EG_Q / 2) 171 sbad[1] = svals[1] 172 gv_check("neg-control-a-CORRUPTED-singular-value-drives-the-SVD-residual-far-above-the-bound" as *u8, 173 (nx_matrix_svd_residual(sa, su, sbad, svv) > EG_RESID_MAX * 4) as i64, ctr) 174 175 // RANK DEFICIENCY IS REPORTED, NOT PAPERED OVER: [[1,1],[1,1]] has rank 1, so the second singular 176 // value must vanish. That zero is how a caller learns the matrix has a null space. 177 let rd: *Matrix = nx_matrix_alloc(2, 2) 178 nx_matrix_set(rd, 0, 0, 1 * EG_Q); nx_matrix_set(rd, 0, 1, 1 * EG_Q) 179 nx_matrix_set(rd, 1, 0, 1 * EG_Q); nx_matrix_set(rd, 1, 1, 1 * EG_Q) 180 let ru: *Matrix = nx_matrix_alloc(2, 2) 181 let rv: *Matrix = nx_matrix_alloc(2, 2) 182 let rvals: *i64 = (sys_mmap(16)) as *i64 183 gv_check_eq("a-RANK-DEFICIENT-matrix-is-accepted" as *u8, 184 nx_matrix_svd(rd, ru, rvals, rv, 12), NXM_SVD_OK, ctr) 185 gv_check_near("a-RANK-DEFICIENT-matrix-returns-a-ZERO-second-singular-value-so-the-null-space-is-visible" as *u8, 186 rvals[1], 0, EG_TOL, ctr) 187 188 let badu: *Matrix = nx_matrix_alloc(3, 3) 189 gv_check_eq("neg-control-NON-CONFORMING-output-dimensions-are-REFUSED-by-their-own-named-rule" as *u8, 190 nx_matrix_svd(sa, badu, svals, svv, 12), NXM_SVD_ERR_DIMS, ctr) 191 192 // ---- EMITTED VALUES ---------------------------------------------------------------------------- 193 gv_values_head() 194 gv_kv("svd_largest_singular_value" as *u8, svals[0]) 195 gv_kv("svd_smallest_singular_value" as *u8, svals[1]) 196 gv_kv("svd_frobenius_from_singular_values" as *u8, frob_sv) 197 gv_kv("svd_frobenius_from_entries" as *u8, frob_a) 198 gv_kv("svd_reconstruction_residual_q14" as *u8, sresid) 199 gv_kv("svd_rank_deficient_second_value" as *u8, rvals[1]) 200 gv_kv("fixed_point_unit_q14" as *u8, EG_Q) 201 gv_kv("known_2x2_larger_eigenvalue" as *u8, hi) 202 gv_kv("known_2x2_smaller_eigenvalue" as *u8, lo) 203 gv_kv("known_2x2_expected_larger" as *u8, 3 * EG_Q) 204 gv_kv("known_2x2_expected_smaller" as *u8, 1 * EG_Q) 205 gv_kv("reconstruction_residual_q14" as *u8, resid) 206 gv_kv("residual_bound_declared" as *u8, EG_RESID_MAX) 207 gv_kv("residual_with_one_corrupted_eigenvalue" as *u8, resid_bad) 208 gv_kv("eigenvector_norm_q14" as *u8, col0) 209 gv_kv("eigenvector_dot_product_q14" as *u8, dot) 210 gv_kv("diagonal_3x3_spectrum_sum" as *u8, evd[0] + evd[1] + evd[2]) 211 gv_kv("singular_matrix_small_eigenvalue" as *u8, small) 212 213 let rcx: i64 = gv_verdict("matrix-eig-gate" as *u8, ctr, 214 // NOTE DELIBERATELY CARRIES NO COUNTS. A verdict note that recites how many teeth or how many 215 // refusals it has is a duplicate ruler beside the tooth list, and it drifts toward understatement 216 // the moment a tooth is added -- the direction nobody audits. Each tooth states itself. 217 "symmetric eigenvalues by cyclic Jacobi and a singular value decomposition built on them, in Q14 fixed point, proved by invariants that hold for any correct implementation -- trace, orthonormality, the Frobenius and determinant identities -- with published reconstruction residuals whose own sensitivity is bite-proven, and every refusal naming the rule that fired" as *u8) 218 sys_exit(rcx) 219 return rcx 220}