We study iterative methods based on Krylov subspaces for low-rank approximation under any Schatten-norm. Here, given access to a matrix A through matrix-vector products, an accuracy parameter , and a target rank , the goal is to find a rank-matrix Z with orthonormal columns such that, where M denotes the ℓ norm of the the singular values of M. For the special cases of = 2 (Frobenius norm) and = ∞ (Spectral norm), Musco and Musco (NeurIPS 2015) obtained an algorithm based on Krylov methods that uses ˜ ( / √ ) matrix-vector products, improving on the naïve ˜ ( / ) dependence obtainable by the power method, where ˜ (•) suppresses poly(log( / )) factors.Our main result is an algorithm that uses only ˜ ( 1/6 / 1/3 ) matrix-vector products, and works for all, not necessarily constant, ≥ 1. For = 2 our bound improves the previous ˜ ( / 1/2 ) bound to ˜ ( / 1/3 ). Since the Schatten-and Schatten-∞ norms of any matrix are the same up to a 1 + factor when ≥ (log )/ , our bound recovers the result of Musco and Musco for = ∞. Further, we prove a matrix-vector query lower bound of Ω(1/ 1/3 ) for any fixed constant ≥ 1, showing that surprisingly Θ(1/ 1/3 ) is the optimal complexity for constant .To obtain our results, we introduce several new techniques, including optimizing over multiple Krylov subspaces simultaneously, and pinching inequalities for partitioned operators. Our lower bound for ∈ [1, 2] uses the Araki-Lieb-Thirring trace inequality, whereas for > 2, we appeal to a norm-compression inequality for aligned partitioned operators. As our algorithms only require matrix-vector product access, they can be applied in settings where alternative techniques such as sketching cannot, e.g., to covariance matrices, Hessians defined implicitly by a neural network, and arbitrary polynomials of a matrix.