Differentiating code

One step in the spectral upsampling method described by Jakob and Hanika (2019) is computing an error function, which is used to find the model parameters. I wanted to provide derivatives of the error function with respect to its arguments to the optimization algorithm I used, in hopes for faster and more reliable convergence.

How does one differentiate Python code though? Every program, no matter how complicated, is a series of simple operations and the chain rule can be used. Computer algebra systems (the one I use is http://maxima.sourceforge.net/Maxima) can help with the process, making this a trivial task.

The error function

The algorithm is as follows:

• Given a triplet of coefficients (model parameters), , , , compute a reflectance according to the model.

• Compute the radiosity of a surface with the reflectance, lit by the given illuminant. This is done by simply multiplying the reflectance by the illuminant's spectral distribution.

• Compute the CIE 1931 XYZ tristimulus values describing how the reflected light would be perceived.

where , and are the cmfs (color matching functions), quantizing human color vision. They are closely related to (but not the same as) spectral sensivities of the three different kinds of cone cells in a human eye. In particular is the luminosity function used in photometry, describing how bright light of different wavelengths appears. The normalization constant ,

assures that for perfectly white surfaces and that arbitrary units can be used for the illuminant and the cmfs.

• Compute the CIE 1976 L*a*b* colorspace coordinates. Start by computing the values of the intermediate lightness function for the three XYZ components.

Likewise for the remaining two components.

• Compute the Euclidean distance between the coordinates above and the target color (given in the L*a*b* colorspace). This is the CIE 1976 color difference (delta E) formula.

Discretization

Because spectral data is typically discrete, that is measured at some chosen wavelengths (e.g. every 10 nanometers in the visible range), the algorithm above operates on arrays rather than continuous functions. Integrals are approximated simply as the corresponding finite sum:

This method turned out to be both fast and accurate enough in this application.