Computer-automated design and discovery have led to high-performance nanophotonic devices with diverse functionalities. However, massively multichannel systems such as metasurfaces controlling many incident angles and photoniccircuit components coupling many waveguide modes still present a challenge. Conventional methods require M in forward simulations and M in adjoint simulations� 2M in simulations in total�to compute the objective function and its gradient for a design involving the response to M in input channels. Here, we develop a formalism that uses the recently proposed augmented partial factorization method to obtain both the objective function and its gradient for a massively multichannel system in a single or a few simulations, achieving over 2 orders of magnitude speedup and reduced memory usage. We use this method to inverse design a metasurface beam splitter that separates the incident light to the target diffraction orders for all incident angles of interest, a key component of the dot projector for 3D sensing. This formalism enables efficient inverse design for a wide range of multichannel optical systems.