Python and Viv · Numerics
Power iteration
What to notice
w @ w(a dot product of two vectors) issum(w * w)in Viv;@there is the matrix product, and a vector times a vector is elementwise.np.ones(n)iszeros(n) + 1.0.np.linalg.eigh(A)returns two arrays; Viv'seig_sym(A, vals, vecs)fills two arrays the caller made, with the eigenvalues ascending aseighgives them. It uses cyclic Jacobi, not LAPACK.vals[-1]," ".join(f"{x:.6f}" for x in v)and the if-expression in the last print are the same in both.- The last bits of the eigenvalue depend on the order of the sums (NumPy calls BLAS), so both programs print ten decimals and whether the two methods agree to 10−12, not the size of the difference.
Python power_iteration.py · 33 lines
Viv power_iteration.viv · 42 lines
py2viv refuses power_iteration.py
1 imports
4import numpy as np
2# nothing to import: @ and eig_sym are built in
py2viv refuses this program:
power_iteration.py:27:14: not numeric Python: np.linalg.eigh's eigenvectors
vals, vecs = np.linalg.eigh(A)
^^^^^^^^^^^^^^^^^
Each eigenvector's sign (and a degenerate pair's basis) is the solver's choice, so Viv's Jacobi and LAPACK disagree; use np.linalg.eigvalsh for the eigenvalues.
power_iteration.py:28:56: not numeric Python: iterating over a value of unknown type
print("all eigenvalues", " ".join(f"{x:.10f}" for x in vals))
^^^^
Fix the refusal that produced it.
2 constructs outside the numeric subset (tools/py2viv/SUBSET.md) 2 a symmetric matrix
7A = np.array([[4.0, 1.0, 0.0],
8 [1.0, 3.0, 1.0],
9 [0.0, 1.0, 2.0]])
5A = [[4.0, 1.0, 0.0],
6 [1.0, 3.0, 1.0],
7 [0.0, 1.0, 2.0]]
(refused)
3 power iteration with the Rayleigh quotient
12def power_iteration(A, iters):
13 v = np.ones(A.shape[0])
14 lam = 0.0
15 for _ in range(iters):
16 w = A @ v
17 v = w / np.sqrt(w @ w)
18 lam = v @ A @ v
19 return lam, v
20
21ITERS = 60
22lam, v = power_iteration(A, ITERS)
23print(f"largest eigenvalue {lam:.10f}")
24print("eigenvector", " ".join(f"{x:.6f}" for x in v))
10def power_iteration(A, iters) {
11 n = shape(A)[0]
12 v = zeros(n) + 1.0
13 lam = 0.0
14 for k to iters {
15 w = A @ v
16 v = w / sqrt(sum(w * w))
17 lam = sum(v * (A @ v))
18 }
19 return lam, v
20}
21
22ITERS = 60
23lam, v = power_iteration(A, ITERS)
24print(f"largest eigenvalue {lam:.10f}")
25print("eigenvector", " ".join(f"{x:.6f}" for x in v))
(refused)
4 against the library eigensolver
27vals, vecs = np.linalg.eigh(A)
28print("all eigenvalues", " ".join(f"{x:.10f}" for x in vals))
29AGREE = 1e-12
30print("agrees with eigh to 1e-12:", "yes" if abs(lam - vals[-1]) < AGREE else "no")
28n = shape(A)[0]
29vals = zeros(n)
30vecs = zeros(n, n)
31eig_sym(A, vals, vecs)
32print("all eigenvalues", " ".join(f"{x:.10f}" for x in vals))
33AGREE = 1e-12
34print("agrees with eigh to 1e-12:", "yes" if abs(lam - vals[-1]) < AGREE else "no")
(refused)
5 checks
33assert np.allclose(A @ v, lam * v, atol=1e-8)
37test "eigenpair" {
38 assert_allclose(A @ v, lam * v, 1e-8)
39}
(refused)
What they print
identical — the two programs print the same text.
Dimmed lines are left out of the comparison. Each Viv test block prints one PASS line; the Python checks the same facts with assert, which prints nothing.
$ python3 power_iteration.py
largest eigenvalue 4.7320508076eigenvector 0.788675 0.577350 0.211325all eigenvalues 1.2679491924 3.0000000000 4.7320508076agrees with eigh to 1e-12: yes
$ viv run power_iteration.viv
largest eigenvalue 4.7320508076eigenvector 0.788675 0.577350 0.211325all eigenvalues 1.2679491924 3.0000000000 4.7320508076agrees with eigh to 1e-12: yesPASS eigenpair