Bundle adjustment is a fundamental problem in computer vision, with important applications such as 3D structure reconstruction from 2D images. This paper focuses on large-scale bundle adjustment tasks, e.g., city-wide 3D reconstruction, which require highly efficient solutions. For this purpose, it is common to apply the Levenberg-Marquardt algorithm, whose bottleneck lies in solving normal equations. The majority of recent methods focus on achieving scalability through modern hardware such as GPUs and distributed systems. On the other hand, the core of the solution, i.e., the math underlying the optimizer for the normal equations, remains largely unimproved since the proposal of the classic parallel bundle adjustment (PBA) algorithm, which increasingly becomes a major limiting factor for the scalability of bundle adjustment. This paper proposes parallel preconditioned conjugate gradient (PPCG) method, a novel parallel method for bundle adjustment based on preconditioned conjugate gradient, which achieves significantly higher efficiency and scalability than existing methods on the algorithmic level. The main idea is to exploit the sparsity of the Hessian matrix and reduce its structure parameters through an effective parallel Schur complement method; the result of this step is then fed into our carefully designed PPCG method which reduces matrix operations that are either expensive (e.g., large matrix reverse or multiplications) or scales poorly to multi-processors (e.g., parallel Reduce operators). Extensive experiments demonstrate that PPCG outperforms existing optimizers by large margins, on a wide range of datasets.