I'm trying to implement the aforementioned procedure using this website as a reference. At the end of the page the algorithm is described as follows:
I think I've mapped the given algorithm to code correctly. However when I use the U and V obtained in the procedure to compute $U^{T}AV$ I don't get a bidiagonal matrix. The values along the diagonals indeed are the ones computed within the procedure but the rest of the values are not zeros. Here's my code:
def golub_kahan(a):
n = a.shape[1]
v = np.ones(n, dtype="float32") / np.sqrt(n)
u = np.zeros(a.shape[0], dtype="float32")
beta = 0
U, V = np.zeros_like(a, dtype="float32"), np.zeros((n,n), dtype="float32")
for i in range(n):
V[:, i] = v
u = a @ v - beta * u
alpha = np.linalg.norm(u)
u /= alpha
U[:, i] = u
v = a.T @ u - alpha * v
beta = np.linalg.norm(v)
v /= beta
return U, V
I'm not sure where my mistake is, did I misunderstand the algorithm somehow?
1e-16
times as large as the other values. $\endgroup$ – wim Dec 28 '18 at 10:16