We present an effective warm-starting scheme for solving large linear complementarity problems (LCPs) arising from Nash equilibrium problems. The approach generates high-quality starting points that, when passed to the PATH solver, yield substantial reductions in computational time and variance. Our warm-start routine reformulates each agent’s LP using strong duality, leading to a master problem with bilinear constraints that is equivalent to the original LCP. Bilinear terms emerge, for example, from system-level variables in the overall equilibrium problem that are exogenous to individual agent optimization problems. Approximate solutions are obtained using the difference-of-convex function algorithm for bilinear terms (DCA-BL) or a spatial branch-and-bound method (SBB). Unlike conventional bilinear approximation schemes, such as McCormick envelopes, DCA-BL does not rely on tight variable bounds. We test the scheme on a realistic LCP instance derived from a stochastic natural gas equilibrium model with nearly 100,000 variables. Without warm starts, PATH struggles to solve these instances within 24 hours. With DCA-BL or SBB warm starts, solution times drop significantly; the largest instance is solved in about one hour after two hours of warm start. While both warm-start approaches yield faster and less variable computational times, experiments suggest that DCA-BL provides the best starting point, as measured by the resulting PATH runtime. Although demonstrated on a specific LCP, the warm-start method extends to any LCP derived from the KKT conditions of LPs for each agent combined with linear system-level constraints.