Surface parameterization is a fundamental problem in computer graphics and geometric processing; thus, it has been widely used in many applications, such as texture mapping, remeshing, intersurface mapping, and shape analysis (cf. extensive surveys [Floater and Hormann 2005; Hormann et al. 2007; Sheffer et al. 2006]). These tasks rely on the computation of a foldover-free, low-distortion parameterization. Linear methods, such as Tutte’s method and its variants [Aigerman and Lipman 2015; Floater 2003; Tutte 1963], provide injectivity-guaranteed parameterization, but usually exhibit extremely large distortion for complex inputs. Nonlinear methods formulate parameterization as an optimization problem by minimizing an energy functions with some constraints that preserve the orientations of triangles. Generally, the objective function includes a low-distortion term that is large when the triangle is highly distorted. It may also include a foldover-preventing term that goes to infinity when a triangle flips or degenerates. These objective functions are highly non-convex and non-linear; thus, they are numerically difficult to optimize, particularly for large-scale inputs.