  We study iterative methods based on Krylov subspaces for low-rank
approximation under any Schatten-$p$ norm. Here, given access to a matrix $A$
through matrix-vector products, an accuracy parameter $\epsilon$, and a target
rank $k$, the goal is to find a rank-$k$ matrix $Z$ with orthonormal columns
such that $\| A(I -ZZ^\top)\|_{S_p} \leq (1+\epsilon)\min_{U^\top U = I_k}
\|A(I - U U^\top)\|_{S_p}$, where $\|M\|_{S_p}$ denotes the $\ell_p$ norm of
the the singular values of $M$. For the special cases of $p=2$ (Frobenius norm)
and $p = \infty$ (Spectral norm), Musco and Musco (NeurIPS 2015) obtained an
algorithm based on Krylov methods that uses $\tilde{O}(k/\sqrt{\epsilon})$
matrix-vector products, improving on the na\"ive $\tilde{O}(k/\epsilon)$
dependence obtainable by the power method, where $\tilde{O}$ suppresses
poly$(\log(dk/\epsilon))$ factors.
  Our main result is an algorithm that uses only
$\tilde{O}(kp^{1/6}/\epsilon^{1/3})$ matrix-vector products, and works for all
$p \geq 1$. For $p = 2$ our bound improves the previous
$\tilde{O}(k/\epsilon^{1/2})$ bound to $\tilde{O}(k/\epsilon^{1/3})$. Since the
Schatten-$p$ and Schatten-$\infty$ norms are the same up to a $(1+
\epsilon)$-factor when $p \geq (\log d)/\epsilon$, our bound recovers the
result of Musco and Musco for $p = \infty$. Further, we prove a matrix-vector
query lower bound of $\Omega(1/\epsilon^{1/3})$ for any fixed constant $p \geq
1$, showing that surprisingly $\tilde{\Theta}(1/\epsilon^{1/3})$ is the optimal
complexity for constant~$k$.
  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 $p \in [1,2]$ uses
the Araki-Lieb-Thirring trace inequality, whereas for $p>2$, we appeal to a
norm-compression inequality for aligned partitioned operators.
