Rank-k reconstruction

Unfold, take the uncentered SVD, truncate, and fold back:

M = tensor.reshape(H, W*C)
U, s, Vt = np.linalg.svd(M, full_matrices=False)
M_k = (U[:, :k] * s[:k]) @ Vt[:k, :]       # broadcasting; no diag matrix needed
restored = M_k.reshape(H, W, C)

U[:, :k] is a row basis; s[:k, None] * Vt[:k] gives k latent scores per column/channel pair.

The result has matrix rank at most k in this unfolding — not necessarily CP rank k. Rank 0 returns zeros, and integer images are promoted to float64 so nothing rounds back to uint8.