As I (Sogo) understand Masayuki's suggestion:
 
Let s_α(λ) the sky spectrum from the α-th fiber.
We first decompose s_α(λ) = sd_α(λ) + sc_α(λ),
where sd_α(λ) is a discrete spectrum and sc_α(λ) is a continuum.
 
Somehow we model sc_α(λ).
 
The discrete spectrum sd_α(λ) will be reconstructible as
    sd_α(λ) = Σ_i a_i P(λ - λ_i; θ_{αi})
if we have a perfect dictionary of emission lines

).
 
After we represent the sky spectrum as
    s_α(λ) = sc_α(λ) + Σ_i a_i P(λ - λ_i; θ_{αi}),
we then somehow tweak sc_α(λ), a_i and θ_{αi} before subtracting
s_α(λ) from an object spectrum, in a way that the sky spectrum
will be removed neatly. This subtraction also includes a
subtle transformation λ ↦ λ', which is almost identity, in order
to absorb errors in wavelength calibration.

edit (2018-11-09)

To be more correct:

0. We assume that all fibers in an exposure share a single, common sky spectrum.

1. From arc exposure(s) we construct a 1-D PSF model for each fiber f:

PSF = PSF(λ - λ_p; f, λ_p).

The PSF depends on the fiber_id f and the peak position λ_p.

We save this PSF model as one of calibration data, just as they save flat and fibertrace.

2. To process spectra (mixture of sky spectra and object spectra),
we first load the PSF model that has been constructed beforehand.
Then, assuming a sky spectrum of fiber f can be expressed as a sum of PSFs:

sky_f(λ) ~ skymodel(λ; (A_p)) = Σ_p A_p PSF(λ - λ_p; f, λ_p),

we determine A_p by minimizing the difference between lhs and rhs.
Note the amplitude A_p depends only on the peak position p.
It does not depend on the fiber f, in accordance with the assumption 0.

3. We subtract the same skymodel(λ; (A_p)) from every object spectra.

Note:

• In 1., we need arc spectra with continua subtracted.

• In 2., we have to separate continua from sky spectra.
  If an upstream module can do this separation, we will just use it.
  Otherwise we have to do it for ourselves.