Performance of the SUGAR model makes it an excellent candidate for use with surveys such as ZTF, LSST or WFIRST, and offers an alternative way of going beyond stretch and color to measure distance with SNIa.

**Acknowledgements.** We thank the technical staff of the University of Hawaii 2.2-m telescope, and Dan Birchall for observing assistance. We recognize the significant cultural role of Mauna Kea within the indigenous Hawaiian community, and we appreciate the opportunity to conduct observations from this revered site. This work was supported in part by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC025CH11231. Support in France was provided by CNRS/IN2P3, CNRS/INSU, and PNC; LPNHE acknowledges support from LABEX ILP, supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11- IDEX-0004-02. NC is grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the University de Lyon for its financial support within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) of the French government operated by the National Research Agency (ANR). Support in Germany was provided by DFG through TRR33 "The Dark Universe" and by DLR through grants FKZ 50OR1503 and FKZ 50OR1602. In China support was provided by Tsinghua University 985 grant and NSFC grant No. 11173017. Some results were obtained using resources and support from the National Energy Research Scientific Computing Center, supported by the Director, Office of Science, Office of Advanced Scientific Computing Research of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank the Gordon and Betty Moore Foundation for their continuing support. Additional support was provided by NASA under the Astrophysics Data Analysis Program grant 15-ADAP15-0256 (PI:Aldering). We also thank the High Performance Research and Education Network (HPWREN), supported by National Science Foundation Grant Nos. 0087344 and 0426879. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 759194 – USNAC). PFL acknowledges support from the National Science Foundation grant PHY-1404070. The work of MVP (participation in SUGAR implementation in sncosmo) was supported by Russian Science Foundation grant 18-72-00159. We thank Claire-Alice Hébert for reviewing and giving helpful advice on this paper.

## References

Aldering, G., Adam, G., Antilogus, P., et al. 2002, SPIE Conf. Ser., 4836, 61

Aldering, G., Antilogus, P., Bailey, S., et al. 2006, ApJ, 650, 510

Amanullah, R., Johansson, J., Goobar, A., et al. 2015, MNRAS, 453, 3300

Arsenijevic, V., Fabbro, S., Mourão, A. M., & Rica da Silva, A. J. 2008, A&A, 492, 535

Astier, P., Guy, J., Regnault, N., et al. 2006, A&A, 447, 31

Bailey, S., Aldering, G., Antilogus, P., et al. 2009, A&A, 500, L17

Bellm, E. 2014, in The Third Hot-wiring the Transient Universe Workshop, eds. P. R. Wozniak, M. J. Graham, A. A. Mahabal, & R. Seaman, 27

Betoule, M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A32

Blondin, S., & Tonry, J. L. 2007, ApJ, 666, 1024

Blondin, S., Mandel, K. S., & Kirshner, R. P. 2011, A&A, 526, A81

Bongard, S., Soulez, F., Thiébaut, É., & Pécontal, É. 2011, MNRAS, 418, 258

Branch, D., Dang, L. C., Hall, N., et al. 2006, PASP, 118, 560

Buton, C., Copin, Y., Aldering, G., et al. 2013, A&A, 549, A8

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Chotard, N. 2011, PhD Thesis, Université Claude Bernard - Lyon I

Chotard, N., Gangler, E., Aldering, G., et al. 2011, A&A, 529, L4

Copeland, E. J., Sami, M., & Tsujikawa, S. 2006, Int. J. Mod. Phys. D, 15, 1753

Delubac, T., Bautista, J. E., Busca, N. G., et al. 2015, A&A, 574, A59

Dempster, A. P., Laird, N. M., & Rubin, D. B. 1977, J. R. Stat. Soc. Ser. B (Methodological), 39, 1

Fakhouri, H. K., Boone, K., Aldering, G., et al. 2015, ApJ, 815, 58

Filippenko, A. V., Richmond, M. W., Matheson, T., et al. 1992a, ApJ, 384, L15

Filippenko, A. V., Richmond, M. W., Branch, D., et al. 1992b, AJ, 104, 1543

Folatelli, G. 2004, New Astron. Rev., 48, 623

Foley, R. J., & Kasen, D. 2011, ApJ, 729, 55

Garavini, G., Folatelli, G., Nobili, S., et al. 2007, A&A, 470, 411

Ghahramani, Z., & Hinton, G. E. 1997, The EM Algorithm for Mixtures of Factor Analyzers, Tech. rep.

Guillochon, J., Parrent, J., Kelley, L. Z., & Margutti, R. 2017, ApJ, 835, 64

Guy, J., Astier, P., Baumont, S., et al. 2007, A&A, 466, 11

Guy, J., Sullivan, M., Conley, A., et al. 2010, A&A, 523, A7

Hamuy, M., Phillips, M. M., Maza, J., et al. 1995, AJ, 109, 1

Kelly, P. L., Hicken, M., Burke, D. L., Mandel, K. S., & Kirshner, R. P. 2010, ApJ, 715, 743

Kim, A. G., Thomas, R. C., Aldering, G., et al. 2013, ApJ, 766, 84

Lantz, B., Aldering, G., Antilogus, P., et al. 2004, SPIE Conf. Ser., 5249, 146

Léget, P. F. 2016, PhD Thesis, Université Blaise Pascal

LSST Dark Energy Science Collaboration 2012, ArXiv e-prints [arXiv:1211.0310]

Maguire, K., Sullivan, M., Ellis, R. S., et al. 2012, MNRAS, 426, 2359

Mandel, K. S., Foley, R. J., & Kirshner, R. P. 2014, ApJ, 797, 75

Mandel, K. S., Scolnic, D. M., Shariff, H., Foley, R. J., & Kirshner, R. P. 2017, ApJ, 842, 93

Nordin, J., Östman, L., Goobar, A., et al. 2011, A&A, 526, A119

Nordin, J., Aldering, G., Antilogus, P., et al. 2018, A&A, 614, A71

Nugent, P., Phillips, M., Baron, E., Branch, D., & Hauschildt, P. 1995, ApJ, 455, L147

Pearson, K. 1901, Phil. Mag., 2, 559

Pereira, R., Thomas, R. C., Aldering, G., et al. 2013, A&A, 554, A27

Perlmutter, S., Gabi, S., Goldhaber, G., et al. 1997, ApJ, 483, 565

Perlmutter, S., Aldering, G., della Valle, M., et al. 1998, Nature, 391, 51

Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565

Phillips, M. M. 1993, ApJ, 413, L105

Planck Collaboration XIII. 2016, A&A, 594, A13

Polin, A., Nugent, P., & Kasen, D. 2019, ApJ, 873, 1

Pskovskii, I. P. 1977, Sov. Astron., 21, 675

Pskovskii, Y. P. 1984, Sov. Astron., 28, 658

Rasmussen, C. E., & Williams, C. K. 2006, Gaussian Processes for Machine Learning (Cambridge: The MIT Press), 38, 715

Rest, A., Scolnic, D., Foley, R. J., et al. 2014, ApJ, 795, 44

Riess, A., Press, W., & Kirshner, R. 1996, AJ, 473, 88

Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009

Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56

Rigault, M., Copin, Y., Aldering, G., et al. 2013, A&A, 560, A66

Rigault, M., Aldering, G., Kowalski, M., et al. 2015, ApJ, 802, 20

Rigault, M., Brinnel, V., Aldering, G., et al. 2018, A&A, submitted

Roman, M., Hardin, D., Betoule, M., et al. 2018, A&A, 615, A68

Rust, B. W. 1974, PhD Thesis, Oak Ridge National Lab., TN

Sasdelli, M., Hillebrandt, W., Aldering, G., et al. 2015, MNRAS, 447, 1247

Saunders, C., Aldering, G., Antilogus, P., et al. 2018, ApJ, 869, 167

Savitzky, A., & Golay, M. J. E. 1964, Anal. Chem., 36, 1627

Scalzo, R. A., Aldering, G., Antilogus, P., et al. 2010, ApJ, 713, 1073

Schmidt, B. P., Suntzeff, N. B., Phillips, M. M., et al. 1998, ApJ, 507, 46

Scolnic, D. M., Riess, A. G., Foley, R. J., et al. 2014, ApJ, 780, 37

Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101

Silverman, J. M., Kong, J. J., & Filippenko, A. V. 2012, MNRAS, 425, 1819

Spearman, C. 1904, Am. J. Psychol., 15, 201

Spearman, C. 1927, The Abilities of Man (New York: Macmillan)

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv e-prints [arXiv:1503.03757]

Sullivan, M., Conley, A., Howell, D. A., et al. 2010, MNRAS, 406, 782

Suzuki, N., Rubin, D., Lidman, C., et al. 2012, ApJ, 746, 85

Tripp, R. 1998, A&A, 331, 815

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528

Wang, X., Filippenko, A. V., Ganeshalingam, M., et al. 2009, ApJ, 699, L139
Performance of the SUGAR model makes it an excellent candidate for use with surveys such as ZTF, LSST or WFIRST, and offers an alternative way of going beyond stretch and color to measure distance with SNIa.

**Acknowledgements.** We thank the technical staff of the University of Hawaii 2.2-m telescope, and Dan Birchall for observing assistance. We recognize the significant cultural role of Mauna Kea within the indigenous Hawaiian community, and we appreciate the opportunity to conduct observations from this revered site. This work was supported in part by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC025CH11231. Support in France was provided by CNRS/IN2P3, CNRS/INSU, and PNC; LPNHE acknowledges support from LABEX ILP, supported by French state funds managed by the ANR within the Investissements d'Avenir programme under reference ANR-11- IDEX-0004-02. NC is grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the University de Lyon for its financial support within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) of the French government operated by the National Research Agency (ANR). Support in Germany was provided by DFG through TRR33 "The Dark Universe" and by DLR through grants FKZ 50OR1503 and FKZ 50OR1602. In China support was provided by Tsinghua University 985 grant and NSFC grant No. 11173017. Some results were obtained using resources and support from the National Energy Research Scientific Computing Center, supported by the Director, Office of Science, Office of Advanced Scientific Computing Research of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank the Gordon and Betty Moore Foundation for their continuing support. Additional support was provided by NASA under the Astrophysics Data Analysis Program grant 15-ADAP15-0256 (PI:Aldering). We also thank the High Performance Research and Education Network (HPWREN), supported by National Science Foundation Grant Nos. 0087344 and 0426879. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 759194 – USNAC). PFL acknowledges support from the National Science Foundation grant PHY-1404070. The work of MVP (participation in SUGAR implementation in sncosmo) was supported by Russian Science Foundation grant 18-72-00159. We thank Claire-Alice Hébert for reviewing and giving helpful advice on this paper.

## References

Aldering, G., Adam, G., Antilogus, P., et al. 2002, SPIE Conf. Ser., 4836, 61

Aldering, G., Antilogus, P., Bailey, S., et al. 2006, ApJ, 650, 510

Amanullah, R., Johansson, J., Goobar, A., et al. 2015, MNRAS, 453, 3300

Arsenijevic, V., Fabbro, S., Mourão, A. M., & Rica da Silva, A. J. 2008, A&A, 492, 535

Astier, P., Guy, J., Regnault, N., et al. 2006, A&A, 447, 31

Bailey, S., Aldering, G., Antilogus, P., et al. 2009, A&A, 500, L17

Bellm, E. 2014, in The Third Hot-wiring the Transient Universe Workshop, eds. P. R. Wozniak, M. J. Graham, A. A. Mahabal, & R. Seaman, 27

Betoule, M., Kessler, R., Guy, J., et al. 2014, A&A, 568, A32

Blondin, S., & Tonry, J. L. 2007, ApJ, 666, 1024

Blondin, S., Mandel, K. S., & Kirshner, R. P. 2011, A&A, 526, A81

Bongard, S., Soulez, F., Thiébaut, É., & Pecontal, É. 2011, MNRAS, 418, 258

Branch, D., Dang, L. C., Hall, N., et al. 2006, PASP, 118, 560

Buton, C., Copin, Y., Aldering, G., et al. 2013, A&A, 549, A8

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Chotard, N. 2011, PhD Thesis, Université Claude Bernard - Lyon I

Chotard, N., Gangler, E., Aldering, G., et al. 2011, A&A, 529, L4

Copeland, E. J., Sami, M., & Tsujikawa, S. 2006, Int. J. Mod. Phys. D, 15, 1753

Delubac, T., Bautista, J. E., Busca, N. G., et al. 2015, A&A, 574, A59

Dempster, A. P., Laird, N. M., & Rubin, D. B. 1977, J. R. Stat. Soc. Ser. B (Methodological), 39, 1

Fakhouri, H. K., Boone, K., Aldering, G., et al. 2015, ApJ, 815, 58

Filippenko, A. V., Richmond, M. W., Matheson, T., et al. 1992a, ApJ, 384, L15

Filippenko, A. V., Richmond, M. W., Branch, D., et al. 1992b, AJ, 104, 1543

Folatelli, G. 2004, New Astron. Rev., 48, 623

Foley, R. J., & Kasen, D. 2011, ApJ, 729, 55

Garavini, G., Folatelli, G., Nobili, S., et al. 2007, A&A, 470, 411

Ghahramani, Z., & Hinton, G. E. 1997, The EM Algorithm for Mixtures of Factor Analyzers, Tech. rep.

Guillochon, J., Parrent, J., Kelley, L. Z., & Margutti, R. 2017, ApJ, 835, 64

Guy, J., Astier, P., Baumont, S., et al. 2007, A&A, 466, 11

Guy, J., Sullivan, M., Conley, A., et al. 2010, A&A, 523, A7

Hamuy, M., Phillips, M. M., Maza, J., et al. 1995, AJ, 109, 1

Kelly, P. L., Hicken, M., Burke, D. L., Mandel, K. S., & Kirshner, R. P. 2010, ApJ, 715, 743

Kim, A. G., Thomas, R. C., Aldering, G., et al. 2013, ApJ, 766, 84

Lantz, B., Aldering, G., Antilogus, P., et al. 2004, SPIE Conf. Ser., 5249, 146

Léget, P. F. 2016, PhD Thesis, Université Blaise Pascal

LSST Dark Energy Science Collaboration 2012, ArXiv e-prints [arXiv:1211.0310]

Maguire, K., Sullivan, M., Ellis, R. S., et al. 2012, MNRAS, 426, 2359

Mandel, K. S., Foley, R. J., & Kirshner, R. P. 2014, ApJ, 797, 75

Mandel, K. S., Scolnic, D. M., Shariff, H., Foley, R. J., & Kirshner, R. P. 2017, ApJ, 842, 93

Nordin, J., Östman, L., Goobar, A., et al. 2011, A&A, 526, A119

Nordin, J., Aldering, G., Antilogus, P., et al. 2018, A&A, 614, A71

Nugent, P., Phillips, M., Baron, E., Branch, D., & Hauschildt, P. 1995, ApJ, 455, L147

Pearson, K. 1901, Phil. Mag., 2, 559

Pereira, R., Thomas, R. C., Aldering, G., et al. 2013, A&A, 554, A27

Perlmutter, S., Gabi, S., Goldhaber, G., et al. 1997, ApJ, 483, 565

Perlmutter, S., Aldering, G., della Valle, M., et al. 1998, Nature, 391, 51

Perlmutter, S., Aldering, G., Goldhaber, G., et al. 1999, ApJ, 517, 565

Phillips, M. M. 1993, ApJ, 413, L105

Planck Collaboration XIII. 2016, A&A, 594, A13

Polin, A., Nugent, P., & Kasen, D. 2019, ApJ, 873, 1

Pskovskii, I. P. 1977, Sov. Astron., 21, 675

Pskovskii, Y. P. 1984, Sov. Astron., 28, 658

Rasmussen, C. E., & Williams, C. K. 2006, Gaussian Processes for Machine Learning (Cambridge: The MIT Press), 38, 715

Rest, A., Scolnic, D., Foley, R. J., et al. 2014, ApJ, 795, 44

Riess, A., Press, W., & Kirshner, R. 1996, AJ, 473, 88

Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009

Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56

Rigault, M., Copin, Y., Aldering, G., et al. 2013, A&A, 560, A66

Rigault, M., Aldering, G., Kowalski, M., et al. 2015, ApJ, 802, 20

Rigault, M., Brinnel, V., Aldering, G., et al. 2018, A&A, submitted

Roman, M., Hardin, D., Betoule, M., et al. 2018, A&A, 615, A68

Rust, B. W. 1974, PhD Thesis, Oak Ridge National Lab., TN

Sasdelli, M., Hillebrandt, W., Aldering, G., et al. 2015, MNRAS, 447, 1247

Saunders, C., Aldering, G., Antilogus, P., et al. 2018, ApJ, 869, 167

Savitzky, A., & Golay, M. J. E. 1964, Anal. Chem., 36, 1627

Scalzo, R. A., Aldering, G., Antilogus, P., et al. 2010, ApJ, 713, 1073

Schmidt, B. P., Suntzeff, N. B., Phillips, M. M., et al. 1998, ApJ, 507, 46

Scolnic, D. M., Riess, A. G., Foley, R. J., et al. 2014, ApJ, 780, 37

Scolnic, D. M., Jones, D. O., Rest, A., et al. 2018, ApJ, 859, 101

Silverman, J. M., Kong, J. J., & Filippenko, A. V. 2012, MNRAS, 425, 1819

Spearman, C. 1904, Am. J. Psychol., 15, 201

Spearman, C. 1927, The Abilities of Man (New York: Macmillan)

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv e-prints [arXiv:1503.03757]

Sullivan, M., Conley, A., Howell, D. A., et al. 2010, MNRAS, 406, 782

Suzuki, N., Rubin, D., Lidman, C., et al. 2012, ApJ, 746, 85

Tripp, R. 1998, A&A, 331, 815

Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528

Wang, X., Filippenko, A. V., Ganeshalingam, M., et al. 2009, ApJ, 699, L139
# Appendix A: Spectral-indicator measurements

As the spectral indicators are a key ingredient of our statistical analysis, we need a robust and automatic algorithm to derive them from the rest-frame spectra and estimate the associated statistical and systematic uncertainties. The main issue to be addressed is the automatic detection of feature boundaries – usually local extrema – as they shift both along the phase and from one supernova to another. Sometimes, the extremum itself cannot be uniquely defined, for example, when there is a mix of several local extrema at the same position or no significant extremum such as in the blue edge of the $O I \lambda 7773$ feature. The photon noise also limit the accuracy when determining the wavelength of the maximum, and may induce a systematic shift due to the nonlinear process involved in finding an extremum. In this section, we give a brief overview of the general and automatic method developed and implemented by Chotard (2011) to measure these spectral indicators and evaluate their corresponding uncertainties. This method shares some similarities with previous analyses (Folatelli 2004; Garavini et al. 2007; Nordin et al. 2011; Blondin et al. 2011; Silverman et al. 2012): it is based on extrema searches in a fixed spectral domain after smoothing the data. Improving on previous studies, our smoothing procedure is based on an optimal determination of the regularization parameters, and more noteworthy, we developed a thorough determination of the uncertainties based on a Monte Carlo procedure taking as input the variance of our signal.

## A.1. Method

Most of the spectral indicators presented here are defined by the position in wavelength and flux of at least one local minimum or maximum, that is, peaks or troughs of the $SN Ia$ spectral features. In order to compute a precise estimate of these local extrema, a Savitsky-Golay (Savitzky & Golay 1964, SG) smoothing is applied to the original spectra before the rebining at $1500 \text{ km s}^{-1}$ described in Sect. 2. Since the SG optimal window depends both on the underlying spectral shape and on the S/N ratio, which varies across wavelength, an independent smoothing is applied to each of the nine spectral zones of interest defined in Table A.1.

Each extremum is then selected on the smoothed spectrum as the local extremum inside a given wavelength range, which can be found in Tables A.2 and A.3 for the equivalent widths and the feature velocities respectively. These wavelength ranges have been trained on a set of ~50 SNe Ia within a phase range of ±5 days around B-band maximum light so that they match the observed diversity of our spectra. The positions in wavelength and flux of these extrema are then used in the spectral-indicator measurements, either directly, for example, for the velocities, or indirectly, for instance, for the equivalent widths through the definition of their pseudo-continuum.

## A.2. Optimal smoothing

The purpose of smoothing, also known as regularization, is to transform the original noisy data in order to get closer on average to the unknown original spectrum, based on some regularity hypothesis. A parameter describing how smooth the final function is has to be introduced, and has to be estimated, either based on physical consideration or deduced from the data themselves. Here we follow the latter approach and describe an optimal way of setting the smoothing parameter given a class of transformations.

<table><caption>**Table A.1.** Independently smoothed spectral regions, with minimal and maximal boudaries $\lambda_{min}$ and $\lambda_{max}$ (in Å).</caption><thead><tr><th>λ<sub>min</sub>–λ<sub>max</sub></th><th>Elements</th></tr></thead><tbody><tr><td>3450–4070</td><td>Ca II H&K</td></tr><tr><td>3850–4150</td><td>Si II λ4131; Co II</td></tr><tr><td>4000–4610</td><td>Mg II triplet</td></tr><tr><td>4350–5350</td><td>Fe II blend</td></tr><tr><td>5060–5700</td><td>S II W</td></tr><tr><td>5500–6050</td><td>Si II λ5972</td></tr><tr><td>5800–6400</td><td>Si II λ6355</td></tr><tr><td>5500–6400</td><td>Si II λ5972; λ6355</td></tr><tr><td>6500–8800</td><td>O I triplet; Ca II IR</td></tr></tbody></table>

<table><caption>**Table A.2.** Definition of the wavelength regions (in Å) where the extrema are expected to be found.</caption><thead><tr><th>Indicators</th><th>λ<sup>b</sup><sub>min</sub> – λ<sup>b</sup><sub>max</sub></th><th>λ<sup>r</sup><sub>min</sub> – λ<sup>r</sup><sub>max</sub></th></tr></thead><tbody><tr><td>Ca II H&K</td><td>3504–3687</td><td>3830–3990</td></tr><tr><td>Si II λ4131</td><td>3830–3990</td><td>4030–4150</td></tr><tr><td>Mg II</td><td>4030–4150</td><td>4450–4650</td></tr><tr><td>Fe λ4800</td><td>4450–4650</td><td>5050–5285</td></tr><tr><td>S II W</td><td>5050–5285</td><td>5500–5681</td></tr><tr><td>Si II λ5972</td><td>5550–5681</td><td>5850–6015</td></tr><tr><td>Si II λ6355</td><td>5850–6015</td><td>6250–6365</td></tr><tr><td>O I λ7773</td><td>7100–7270</td><td>7720–8000</td></tr><tr><td>Ca II IR</td><td>7720–8000</td><td>8300–8800</td></tr></tbody></table>

**Notes.** The b and r exponents respectively represent the left (blue) and right (red) peak.

<table><caption>**Table A.3.** Wavelength regions (in Å) used to compute the feature velocities, together with their rest-frame wavelengths, $\lambda_0$.</caption><thead><tr><th>Velocity</th><th>λ<sub>min</sub>–λ<sub>max</sub></th><th>λ<sub>0</sub></th></tr></thead><tbody><tr><td>ν(Si II λ4131)</td><td>3963–4034</td><td>4131</td></tr><tr><td>ν(S II W λ5454)</td><td>5200–5350</td><td>5454</td></tr><tr><td>ν(S II W λ5640)</td><td>5351–5550</td><td>5640</td></tr><tr><td>ν(Si II λ6355)</td><td>6000–6210</td><td>6355</td></tr></tbody></table>

### A.2.1. Formalism

Smoothing a noisy spectrum $Y$ consists of the determination of the smoothed spectrum $Y'$, which is a function of the initial flux $Y' = f(Y)$. $Y'$ is an estimator of the noise-free and unknown original flux, convolved by the instrumental resolution: $\hat{Y}$. It is related to the observed flux by:

$$Y = \hat{Y} + N, \qquad (A.1)$$

where $N$ is a realization of the spectral noise vector. The goodness of the smoothing being represented by the error function,

$$\|Y' - \hat{Y}\| \approx (Y' - \hat{Y})^T W (Y' - \hat{Y}), \qquad (A.2)$$

where $W$ is the inverse of the noise covariance matrix, this requires minimization of a quantity which depends both on the noise properties and the signal shape. As the true spectrum, $\hat{Y}$, is unknown, we need to build an estimator of $\|Y' - \hat{Y}\|$ which is independent of $\hat{Y}$. For a linear regularization, one can write:

$$Y' = f(Y) = B_p Y, \qquad (A.3)$$
# Appendix A: Spectral-indicator measurements

As the spectral indicators are a key ingredient of our statistical analysis, we need a robust and automatic algorithm to derive them from the rest-frame spectra and estimate the associated statistical and systematic uncertainties. The main issue to be addressed is the automatic detection of feature boundaries – usually local extrema – as they shift both along the phase and from one supernova to another. Sometimes, the extremum itself cannot be uniquely defined, for example, when there is a mix of several local extrema at the same position or no significant extremum such as in the blue edge of the $O I \lambda 7773$ feature. The photon noise also limit the accuracy when determining the wavelength of the maximum, and may induce a systematic shift due to the nonlinear process involved in finding an extremum. In this section, we give a brief overview of the general and automatic method developed and implemented by Chotard (2011) to measure these spectral indicators and evaluate their corresponding uncertainties. This method shares some similarities with previous analyses (Folatelli 2004; Garavini et al. 2007; Nordin et al. 2011; Blondin et al. 2011; Silverman et al. 2012): it is based on extrema searches in a fixed spectral domain after smoothing the data. Improving on previous studies, our smoothing procedure is based on an optimal determination of the regularization parameters, and more noteworthy, we developed a thorough determination of the uncertainties based on a Monte Carlo procedure taking as input the variance of our signal.

## A.1. Method

Most of the spectral indicators presented here are defined by the position in wavelength and flux of at least one local minimum or maximum, that is, peaks or troughs of the $SN Ia$ spectral features. In order to compute a precise estimate of these local extrema, a Savitsky-Golay (Savitzky & Golay 1964, SG) smoothing is applied to the original spectra before the rebining at $1500 \text{ km s}^{-1}$ described in Sect. 2. Since the SG optimal window depends both on the underlying spectral shape and on the S/N ratio, which varies across wavelength, an independent smoothing is applied to each of the nine spectral zones of interest defined in Table A.1.

Each extremum is then selected on the smoothed spectrum as the local extremum inside a given wavelength range, which can be found in Tables A.2 and A.3 for the equivalent widths and the feature velocities respectively. These wavelength ranges have been trained on a set of ~50 SNe Ia within a phase range of ±5 days around B-band maximum light so that they match the observed diversity of our spectra. The positions in wavelength and flux of these extrema are then used in the spectral-indicator measurements, either directly, for example, for the velocities, or indirectly, for instance, for the equivalent widths through the definition of their pseudo-continuum.

## A.2. Optimal smoothing

The purpose of smoothing, also known as regularization, is to transform the original noisy data in order to get closer on average to the unknown original spectrum, based on some regularity hypothesis. A parameter describing how smooth the final function is has to be introduced, and has to be estimated, either based on physical consideration or deduced from the data themselves. Here we follow the latter approach and describe an optimal way of setting the smoothing parameter given a class of transformations.

<table><caption>**Table A.1.** Independently smoothed spectral regions, with minimal and maximal boudaries $\lambda_{min}$ and $\lambda_{max}$ (in Å).</caption><thead><tr><th>λ<sub>min</sub>–λ<sub>max</sub></th><th>Elements</th></tr></thead><tbody><tr><td>3450–4070</td><td>Ca II H&K</td></tr><tr><td>3850–4150</td><td>Si II λ4131; Co II</td></tr><tr><td>4000–4610</td><td>Mg II triplet</td></tr><tr><td>4350–5350</td><td>Fe II blend</td></tr><tr><td>5060–5700</td><td>S II W</td></tr><tr><td>5500–6050</td><td>Si II λ5972</td></tr><tr><td>5800–6400</td><td>Si II λ6355</td></tr><tr><td>5500–6400</td><td>Si II λ5972; λ6355</td></tr><tr><td>6500–8800</td><td>O I triplet; Ca II IR</td></tr></tbody></table>

<table><caption>**Table A.2.** Definition of the wavelength regions (in Å) where the extrema are expected to be found.</caption><thead><tr><th>Indicators</th><th>λ<sup>b</sup><sub>min</sub>–λ<sup>b</sup><sub>max</sub></th><th>λ<sup>r</sup><sub>min</sub>–λ<sup>r</sup><sub>max</sub></th></tr></thead><tbody><tr><td>Ca II H&K</td><td>3504–3687</td><td>3830–3990</td></tr><tr><td>Si II λ4131</td><td>3830–3990</td><td>4030–4150</td></tr><tr><td>Mg II</td><td>4030–4150</td><td>4450–4650</td></tr><tr><td>Fe λ4800</td><td>4450–4650</td><td>5050–5285</td></tr><tr><td>S II W</td><td>5050–5285</td><td>5500–5681</td></tr><tr><td>Si II λ5972</td><td>5550–5681</td><td>5850–6015</td></tr><tr><td>Si II λ6355</td><td>5850–6015</td><td>6250–6365</td></tr><tr><td>O I λ7773</td><td>7100–7270</td><td>7720–8000</td></tr><tr><td>Ca II IR</td><td>7720–8000</td><td>8300–8800</td></tr></tbody></table>

**Notes.** The *b* and *r* exponents respectively represent the left (blue) and right (red) peak.

<table><caption>**Table A.3.** Wavelength regions (in Å) used to compute the feature velocities, together with their rest-frame wavelengths, $\lambda_0$.</caption><thead><tr><th>Velocity</th><th>λ<sub>min</sub>–λ<sub>max</sub></th><th>λ<sub>0</sub></th></tr></thead><tbody><tr><td>ν(Si II λ4131)</td><td>3963–4034</td><td>4131</td></tr><tr><td>ν(S II W λ5454)</td><td>5200–5350</td><td>5454</td></tr><tr><td>ν(S II W λ5640)</td><td>5351–5550</td><td>5640</td></tr><tr><td>ν(Si II λ6355)</td><td>6000–6210</td><td>6355</td></tr></tbody></table>

### A.2.1. Formalism

Smoothing a noisy spectrum $Y$ consists of the determination of the smoothed spectrum $Y'$, which is a function of the initial flux $Y' = f(Y)$. $Y'$ is an estimator of the noise-free and unknown original flux, convolved by the instrumental resolution: $\hat{Y}$. It is related to the observed flux by:

$$Y = \hat{Y} + N, \qquad (A.1)$$

where $N$ is a realization of the spectral noise vector. The goodness of the smoothing being represented by the error function,

$$||Y' - \hat{Y}|| \approx (Y' - \hat{Y})^T W (Y' - \hat{Y}), \qquad (A.2)$$

where $W$ is the inverse of the noise covariance matrix, this requires minimization of a quantity which depends both on the noise properties and the signal shape. As the true spectrum, $\hat{Y}$, is unknown, we need to build an estimator of $||Y' - \hat{Y}||$ which is independent of $\hat{Y}$. For a linear regularization, one can write:

$$Y' = f(Y) = B_p Y, \qquad (A.3)$$
where $B_p$ is a smoothing matrix which only depends on the chosen smoothing technique and a smoothing parameter $p$. Introducing Eqs. (A.1) and (A.3) into (A.2), one can show that:

$$ \|Y' - \hat{Y}\| = (Y' - Y + N)^T W(Y' - Y + N) \quad (A.4) $$

$$ = \|Y' - Y\| + 2N^T W(B_p - I)(\hat{Y} + N) + N^T W N. \quad (A.5) $$

An estimator of the error function, $\epsilon$, that has to be minimized can be constructed by noticing that

$$ E[\|Y' - \hat{Y}\|] = E[\|Y' - Y\|] + 2 \mathrm{Tr}(B_p) - n, \quad (A.6) $$

where $E$ is the mathematical expectation value, $\mathrm{Tr}$ is the Trace operator, and $n$ is the rank of the vector $Y$. As we have only one realization of $Y$, this translates to:

$$ \epsilon = \|Y' - Y\| + 2 \mathrm{Tr}(B_p) - n, \quad (A.7) $$

which is the quantity that is minimized with respect to $p$.

Given a class of linear smoothing methods $B_p$, we have thus defined a procedure to find the value of $p$ for which the smoothed spectrum best reproduces the original unknown one given an observed spectrum. This method is not exempt from possible overtraining, however this is mitigated by restricting the optimization to a single parameter. This was tested with Monte Carlo simulations.

### A.2.2. B matrix estimation

Among all possibilities, we chose a Savitzky-Golay regularization (Savitzky & Golay 1964, SG) as it allows the reduction of the high frequency noise, while keeping the original shape of the feature. This method relies on fitting a k-order polynomial function for each point $i$ in a fixed window size $p$ (with $p > k + 1$) centered on the current point and is designed to preserve the locus of maxima for even values of $k$. Each observed data point is replaced by its fit value, and the window then moves to the next data point until the spectrum is completely smoothed. As we have fixed the degree of the polynomial function to $k = 2$, the only parameter that has to be optimized is the window size, $p$. This transformation is linear, and the optimal window size $p$ and corresponding smoothing matrix $B_p$ are found by minimizing Eq. (A.7). The value of $p$ is estimated individually for each spectral region of a given spectrum, and the spectral indicator measurements are performed on the corresponding smoothed spectrum, as indicated above.

### A.3. Uncertainties

The uncertainty on a given spectral-indicator measurement arises from the statistical noise of the data and the induced uncertainties on the measurement method parameters. In our approach, these two uncertainties are independently measured using the smoothed spectrum as a reference: the SG window size is typically large enough so that the residual noise can be neglected in the simulations.

A random noise matching the statistical properties of the observed data is then added in order to generate a mock spectrum, and the spectral indicators are measured on this new spectrum as if it were the observed one. After one thousand generations, we are able to derive the statistical fluctuations of the obtained values, which we quote as the statistical uncertainty. This is computed as the standard deviation from the value measured on the real spectrum, thus taking into account a potential

bias. This bias typically corresponds to 10–20% of the total error budget.

In order to save computation time, the initial SG window size $p$ is determined once and for all using the original spectrum. It is then kept constant when measuring spectral indicators on all the simulated spectra. However, the noise affecting the initial spectrum induces an uncertainty in the smoothing parameter. This has been studied on a reduced set of simulations for which the optimal $p$ was derived for each realization of the noise, and the corresponding uncertainty was found to range from 15% to 20% depending on the spectral zone. We then propagate this uncertainty by computing the spectral-indicator values on the original spectrum for several values of $p$ in this 15–20% range. The standard deviation of the resulting spectral-indicator value distribution gives an estimate of the systematic error introduced by the arbitrariness of the smoothing method, which is found to be ~20% of the total error on average. These two uncertainties are quadratically added together to derive the final uncertainty for a given spectral indicator. This estimate takes into account the statistical noise of the spectrum, as well as the induced scatter in wavelength and flux of the extrema. Our measurement errors thus include all nonlinear effects due to limited signal to noise and can be trusted for subsequent statistical analysis. More details could be found in Chotard (2011) and Nordin et al. (2011).

### A.4. Performances

#### A.4.1. Failure rate

When one of the extrema defining a spectral indicator lies on a flat and/or noisy section of the spectrum, its measurement has a chance of failure. In that case, the measurement is automatically rejected and a visual scan using control plots is performed to confirm the actual lack of an extremum. The rejection of the “bad” measurements is performed using a $3\sigma$ clipping in the measured-uncertainty space: if the uncertainty made on a given measurement is larger than $m + 3 \times \text{std}$, where $m$ and $\text{std}$ are the average and standard deviation of this spectral-indicator uncertainty distribution (for the whole sample), the corresponding measurement is rejected. Considering all the spectral indicators measured on the 113 input spectra in a range of phase of $\pm 2.5$ days around maximum light and presented in this paper, the global failure rate of the measurement procedure is less than $10^{-3}$ for the automatic selections mentioned above. If it does happen, the spectral indicators are set to the average value and assigned an infinite error.

#### A.4.2. Quoted uncertainties

A simple test has been performed to confirm the robustness of the method. In our selected sample, 18 SNe Ia have two or more spectra taken in the same night in a phase interval of $\pm 5$ days around maximum light. We then computed the distribution of the pull,

$$ \frac{\delta I}{\sigma_{\delta I}} = \frac{I_1 - I_2}{\sqrt{\sigma_{I_1}^2 + \sigma_{I_2}^2}}, \quad (A.8) $$

for all the spectral indicators $I$ (feature velocity and absorption ratio) measured on each spectrum of a same night (for a same supernova). This distribution is centered around $-0.01$ with a dispersion of $1.06$ which indicates that our estimation of the uncertainty is valid up to a possible underestimation of the error
where $B_p$ is a smoothing matrix which only depends on the chosen smoothing technique and a smoothing parameter $p$. Introducing Eqs. (A.1) and (A.3) into (A.2), one can show that:

$$ \|Y' - \hat{Y}\| = (Y' - Y + N)^T W(Y' - Y + N) \quad (A.4) $$

$$ = \|Y' - Y\| + 2N^T W(B_p - I)(\hat{Y} + N) + N^T W N. \quad (A.5) $$

An estimator of the error function, $\epsilon$, that has to be minimized can be constructed by noticing that

$$ E[\|Y' - \hat{Y}\|] = E[\|Y' - Y\|] + 2 \mathrm{Tr}(B_p) - n, \quad (A.6) $$

where $E$ is the mathematical expectation value, $\mathrm{Tr}$ is the Trace operator, and $n$ is the rank of the vector $Y$. As we have only one realization of $Y$, this translates to:

$$ \epsilon = \|Y' - Y\| + 2 \mathrm{Tr}(B_p) - n, \quad (A.7) $$

which is the quantity that is minimized with respect to $p$.

Given a class of linear smoothing methods $B_p$, we have thus defined a procedure to find the value of $p$ for which the smoothed spectrum best reproduces the original unknown one given an observed spectrum. This method is not exempt from possible overtraining, however this is mitigated by restricting the optimization to a single parameter. This was tested with Monte Carlo simulations.

### A.2.2. B matrix estimation

Among all possibilities, we chose a Savitzky-Golay regularization (Savitzky & Golay 1964, SG) as it allows the reduction of the high frequency noise, while keeping the original shape of the feature. This method relies on fitting a $k$-order polynomial function for each point $i$ in a fixed window size $p$ (with $p > k + 1$) centered on the current point and is designed to preserve the locus of maxima for even values of $k$. Each observed data point is replaced by its fit value, and the window then moves to the next data point until the spectrum is completely smoothed. As we have fixed the degree of the polynomial function to $k = 2$, the only parameter that has to be optimized is the window size, $p$. This transformation is linear, and the optimal window size $p$ and corresponding smoothing matrix $B_p$ are found by minimizing Eq. (A.7). The value of $p$ is estimated individually for each spectral region of a given spectrum, and the spectral indicator measurements are performed on the corresponding smoothed spectrum, as indicated above.

### A.3. Uncertainties

The uncertainty on a given spectral-indicator measurement arises from the statistical noise of the data and the induced uncertainties on the measurement method parameters. In our approach, these two uncertainties are independently measured using the smoothed spectrum as a reference: the SG window size is typically large enough so that the residual noise can be neglected in the simulations.

A random noise matching the statistical properties of the observed data is then added in order to generate a mock spectrum, and the spectral indicators are measured on this new spectrum as if it were the observed one. After one thousand generations, we are able to derive the statistical fluctuations of the obtained values, which we quote as the statistical uncertainty. This is computed as the standard deviation from the value measured on the real spectrum, thus taking into account a potential

bias. This bias typically corresponds to 10–20% of the total error budget.

In order to save computation time, the initial SG window size $p$ is determined once and for all using the original spectrum. It is then kept constant when measuring spectral indicators on all the simulated spectra. However, the noise affecting the initial spectrum induces an uncertainty in the smoothing parameter. This has been studied on a reduced set of simulations for which the optimal $p$ was derived for each realization of the noise, and the corresponding uncertainty was found to range from 15% to 20% depending on the spectral zone. We then propagate this uncertainty by computing the spectral-indicator values on the original spectrum for several values of $p$ in this 15–20% range. The standard deviation of the resulting spectral-indicator value distribution gives an estimate of the systematic error introduced by the arbitrariness of the smoothing method, which is found to be ~20% of the total error on average. These two uncertainties are quadratically added together to derive the final uncertainty for a given spectral indicator. This estimate takes into account the statistical noise of the spectrum, as well as the induced scatter in wavelength and flux of the extrema. Our measurement errors thus include all nonlinear effects due to limited signal to noise and can be trusted for subsequent statistical analysis. More details could be found in Chotard (2011) and Nordin et al. (2011).

### A.4. Performances

#### A.4.1. Failure rate

When one of the extrema defining a spectral indicator lies on a flat and/or noisy section of the spectrum, its measurement has a chance of failure. In that case, the measurement is automatically rejected and a visual scan using control plots is performed to confirm the actual lack of an extremum. The rejection of the “bad” measurements is performed using a $3\sigma$ clipping in the measured-uncertainty space: if the uncertainty made on a given measurement is larger than $m + 3 \times \text{std}$, where $m$ and $\text{std}$ are the average and standard deviation of this spectral-indicator uncertainty distribution (for the whole sample), the corresponding measurement is rejected. Considering all the spectral indicators measured on the 113 input spectra in a range of phase of $\pm 2.5$ days around maximum light and presented in this paper, the global failure rate of the measurement procedure is less than $10^{-3}$ for the automatic selections mentioned above. If it does happen, the spectral indicators are set to the average value and assigned an infinite error.

#### A.4.2. Quoted uncertainties

A simple test has been performed to confirm the robustness of the method. In our selected sample, 18 SNe Ia have two or more spectra taken in the same night in a phase interval of $\pm 5$ days around maximum light. We then computed the distribution of the pull,

$$ \frac{\delta I}{\sigma_{\delta I}} = \frac{I_1 - I_2}{\sqrt{\sigma_{I_1}^2 + \sigma_{I_2}^2}}, \quad (A.8) $$

for all the spectral indicators $I$ (feature velocity and absorption ratio) measured on each spectrum of a same night (for a same supernova). This distribution is centered around $-0.01$ with a dispersion of $1.06$ which indicates that our estimation of the uncertainty is valid up to a possible underestimation of the error
by 6%. This number is small enough so that we can trust our uncertainty estimation for the main analysis.

## Appendix B: Expectation-maximization factor analysis

Dimensionality reduction in the presence of noisy data is often an overlooked problem. However, standard methods like PCA tend to fail at capturing the intrinsic variability of the data and the principal components align with the direction of the noise when the latter becomes important. Factor analysis on the other hand is a statistical method designed to model the covariance structure of high dimensional data using a small number of latent variables. It estimates both the natural variability of the sample and the noise arising from the measurements, under the assumption that the statistics are the same for all data records. Our case is slightly different: on one hand, the noise statistics are different for each measurement, but on the other hand, their variance is already known. We thus adapted the expectation-minimization algorithm presented in Ghahramani & Hinton (1997) to accommodate for the specifics of our problem. The resulting method is also known as probabilistic PCA. The formalism is the following: $x_i$ is a vector of rank $l$ representing the measurement $i$. It is linked to the factor $q_i$, vector of rank $k \le l$ and the noise $\eta_i$ by:

$$x_i = x^0 + \Lambda q_i + \eta_i, \quad (B.1)$$

where $x^0$ is a central value which can be further neglected without loss of generality ($x^0 = 0$), $q_i$ is assumed to follow a normal distribution of unit variance, and $\eta_i$ follows a multivariate normal distribution of variance $\Psi_i$. In our case, $\Psi_i$ is diagonal, a property that can be used to speed-up computations. $\Lambda$ is the matrix containing the $k$ explicative vectors that we need to determine. While $\Lambda$ is not uniquely defined, $\Lambda\Lambda^T$ is and represents the intrinsic covariance of the data, that is, the one we would observe in the absence of noise. The eigenvectors of $\Lambda\Lambda^T$ thus correspond to the $k$ first eigenvectors that PCA would have found in the absence of noise.

To find $\Lambda$, instead of directly maximizing the likelihood of observing $x_i$, the expectation-maximization algorithm introduces the latent variable $q_i$ and then maximizes the expected likelihood over $q_i$. The joint probability of $x_i$ and $q_i$ is the following multivariate normal distribution:

$$P\left(\begin{bmatrix} x_i \\ q_i \end{bmatrix}\right) = N\left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \Lambda\Lambda^T + \Psi_i & \Lambda \\ \Lambda^T & I \end{bmatrix}\right). \quad (B.2)$$

The block-diagonal elements of the covariance matrix represent, respectively, the covariances of $x_i$ and $q_i$, and the non block-diagonal elements represent the covariance arising from the relation (B.1). Expectation-maximization is an iterative procedure which ensures that the likelihood increases at each iteration and it has been shown that the convergence is faster than using a gradient method (Dempster et al. 1977). Each iteration proceeds in two steps. The first step, called the E-step, consists of calculating the expectation of the conditional first and second moments of $q_i$ for a given $\Lambda$:

E-step:

$$q_i \hat{=} E[q|x_i] = \Lambda^T (\Psi_i + \Lambda\Lambda^T)^{-1} x_i \quad (B.3)$$

$$E[q q^T | x_i] = I - \Lambda^T (\Psi_i + \Lambda\Lambda^T)^{-1} \Lambda + q_i q_i^T. \quad (B.4)$$

Once these two quantities are computed, the second step, or M-step, consists of estimating the $\Lambda$ matrix by maximizing the likelihood expectation, which provides the condition:

$$\sum_{i=0}^{N} \Psi_i^{-1} \Lambda E[q q^T | x_i] = \sum_{i=0}^{N} \Psi_i^{-1} x_i q_i^T. \quad (B.5)$$

In the case where $\Psi_i$ is diagonal, one can independently calculate each row of $\Lambda$, noted $\Lambda^j$, according to the relation M-step:

$$\Lambda^j = \left[ \sum_{i=0}^{N} \frac{x_i^j}{\psi_i^{jj}} q_i^T \right] \left[ \sum_{i=0}^{N} \frac{1}{\psi_i^{jj}} E[q q^T | x_i] \right]^{-1}, \quad (B.6)$$

After the last iteration, the internal degeneracies of the description are lifted with the transformation of $\Lambda$ into an orthogonal $\Lambda'$ matrix which satisfies the condition:

$$\Lambda' \Lambda'^T = \Lambda \Lambda^T. \quad (B.7)$$

The columns of $\Lambda'$ are aligned with the eigenvectors of $\Lambda\Lambda^T$ and the square of their norm are the respective eigenvalues.

## Appendix C: Orthogonal distance regression

### C.1. Expectation and maximization steps

Similarly to the factor analysis, $h_i$ and $\Lambda$ are estimated iteratively. The first step consists of minimizing the $\chi^2$ with respect to $h_i$, which amounts to solving the equation:

$$\frac{\partial \chi^2}{\partial h_i} = 0 \quad (C.1)$$

and consequently gives:

E-step:

$$h_i = (\Lambda^T W_{M_i} \Lambda + W_{x_i})^{-1} (\Lambda^T W_{M_i} M_i + W_{x_i} x_i). \quad (C.2)$$

Once the E-step is performed we can estimate the matrix $\Lambda$ by

$$\frac{\partial \chi^2}{\partial \Lambda} = 0, \quad (C.3)$$

which gives:

M-step:

$$\tilde{\Lambda} = \left( \left[ \sum_{i=0}^{N} h_i h_i^T \otimes W_{M_i} \right]^{-1} \left[ \sum_{i=0}^{N} W_{M_i} M_i h_i^T \right] \right)^{-1}, \quad (C.4)$$

where $\tilde{\Lambda}$ is the vector that contains all columns of $\Lambda$ put end-to-end and the $\otimes$ operator is the Kronecker product. The estimation of the $\Lambda$ matrix and the orthogonal projections $h_i$ is done by reiterating the E-step and M-step until $\chi^2$ convergences.

### C.2. Fixing the degeneracies

The free parameters of Eqs. (5) and (22) that we find through orthogonal distance regression contain degenerate degrees of freedom. These must be fixed in order to present an unique an interpretable solution.
by 6%. This number is small enough so that we can trust our uncertainty estimation for the main analysis.

## Appendix B: Expectation-maximization factor analysis

Dimensionality reduction in the presence of noisy data is often an overlooked problem. However, standard methods like PCA tend to fail at capturing the intrinsic variability of the data and the principal components align with the direction of the noise when the latter becomes important. Factor analysis on the other hand is a statistical method designed to model the covariance structure of high dimensional data using a small number of latent variables. It estimates both the natural variability of the sample and the noise arising from the measurements, under the assumption that the statistics are the same for all data records. Our case is slightly different: on one hand, the noise statistics are different for each measurement, but on the other hand, their variance is already known. We thus adapted the expectation-minimization algorithm presented in Ghahramani & Hinton (1997) to accommodate for the specifics of our problem. The resulting method is also known as probabilistic PCA. The formalism is the following: $\mathbf{x}_i$ is a vector of rank $l$ representing the measurement $i$. It is linked to the factor $\mathbf{q}_i$, vector of rank $k \le l$ and the noise $\eta_i$ by:

$$ \mathbf{x}_i = \mathbf{x}^0 + \mathbf{\Lambda}\mathbf{q}_i + \eta_i, \qquad (B.1) $$

where $\mathbf{x}^0$ is a central value which can be further neglected without loss of generality ($\mathbf{x}^0 = 0$), $\mathbf{q}_i$ is assumed to follow a normal distribution of unit variance, and $\eta_i$ follows a multivariate normal distribution of variance $\mathbf{\Lambda}\mathbf{q}_i$. In our case, $\mathbf{\Lambda}$ is diagonal, a property that can be used to speed-up computations. $\mathbf{\Lambda}$ is the matrix containing the $k$ explicative vectors that we need to determine. While $\mathbf{\Lambda}$ is not uniquely defined, $\mathbf{\Lambda}\mathbf{\Lambda}^T$ is and represents the intrinsic covariance of the data, that is, the one we would observe in the absence of noise. The eigenvectors of $\mathbf{\Lambda}\mathbf{\Lambda}^T$ thus correspond to the $k$ first eigenvectors that PCA would have found in the absence of noise.

To find $\mathbf{\Lambda}$, instead of directly maximizing the likelihood of observing $\mathbf{x}_i$, the expectation-maximization algorithm introduces the latent variable $\mathbf{q}_i$ and then maximizes the expected likelihood over $\mathbf{q}_i$. The joint probability of $\mathbf{x}_i$ and $\mathbf{q}_i$ is the following multivariate normal distribution:

$$ P\left(\begin{bmatrix} \mathbf{x}_i \\ \mathbf{q}_i \end{bmatrix}\right) = \mathcal{N}\left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \mathbf{\Lambda}\mathbf{\Lambda}^T + \mathbf{\Psi}_i & \mathbf{\Lambda} \\ \mathbf{\Lambda}^T & \mathbf{I} \end{bmatrix}\right). \qquad (B.2) $$

The block-diagonal elements of the covariance matrix represent, respectively, the covariances of $\mathbf{x}_i$ and $\mathbf{q}_i$, and the non block-diagonal elements represent the covariance arising from the relation (B.1). Expectation-maximization is an iterative procedure which ensures that the likelihood increases at each iteration and it has been shown that the convergence is faster than using a gradient method (Dempster et al. 1977). Each iteration proceeds in two steps. The first step, called the E-step, consists of calculating the expectation of the conditional first and second moments of $\mathbf{q}_i$ for a given $\mathbf{\Lambda}$:

E-step:

$$ \mathbf{q}_i \hat{=} E[\mathbf{q}|\mathbf{x}_i] = \mathbf{\Lambda}^T(\mathbf{\Psi}_i + \mathbf{\Lambda}\mathbf{\Lambda}^T)^{-1}\mathbf{x}_i \qquad (B.3) $$

$$ E[\mathbf{q} \mathbf{q}^T | \mathbf{x}_i] = \mathbf{I} - \mathbf{\Lambda}^T (\mathbf{\Psi}_i + \mathbf{\Lambda}\mathbf{\Lambda}^T)^{-1} \mathbf{\Lambda} + \mathbf{q}_i \mathbf{q}_i^T. \qquad (B.4) $$

Once these two quantities are computed, the second step, or M-step, consists of estimating the $\mathbf{\Lambda}$ matrix by maximizing the likelihood expectation, which provides the condition:

$$ \sum_{i=0}^{N} \Psi_i^{-1} \mathbf{\Lambda} E[\mathbf{q} \mathbf{q}^T | \mathbf{x}_i] = \sum_{i=0}^{N} \Psi_i^{-1} \mathbf{x}_i \mathbf{q}_i^T. \qquad (B.5) $$

In the case where $\mathbf{\Psi}_i$ is diagonal, one can independently calculate each row of $\mathbf{\Lambda}$, noted $\mathbf{\Lambda}^j$, according to the relation M-step:

$$ \mathbf{\Lambda}^j = \left[ \sum_{i=0}^{N} \frac{\mathbf{x}_i^j}{\psi_i^{jj}} \mathbf{q}_i^T \right] \left[ \sum_{i=0}^{N} \frac{1}{\psi_i^{jj}} E[\mathbf{q} \mathbf{q}^T | \mathbf{x}_i] \right]^{-1}, \qquad (B.6) $$

After the last iteration, the internal degeneracies of the description are lifted with the transformation of $\mathbf{\Lambda}$ into an orthogonal $\mathbf{\Lambda}'$ matrix which satisfies the condition:

$$ \mathbf{\Lambda}'\mathbf{\Lambda}'^T = \mathbf{\Lambda}\mathbf{\Lambda}^T. \qquad (B.7) $$

The columns of $\mathbf{\Lambda}'$ are aligned with the eigenvectors of $\mathbf{\Lambda}\mathbf{\Lambda}^T$ and the square of their norm are the respective eigenvalues.

## Appendix C: Orthogonal distance regression

### C.1. Expectation and maximization steps

Similarly to the factor analysis, $\mathbf{h}_i$ and $\mathbf{\Lambda}$ are estimated iteratively. The first step consists of minimizing the $\chi^2$ with respect to $\mathbf{h}_i$, which amounts to solving the equation:

$$ \frac{\partial \chi^2}{\partial \mathbf{h}_i} = 0 \qquad (C.1) $$

and consequently gives:

E-step:

$$ \mathbf{h}_i = (\mathbf{\Lambda}^T \mathbf{W}_{\mathbf{M}_i} \mathbf{\Lambda} + \mathbf{W}_{\mathbf{x}_i})^{-1} (\mathbf{\Lambda}^T \mathbf{W}_{\mathbf{M}_i} \mathbf{M}_i + \mathbf{W}_{\mathbf{x}_i} \mathbf{x}_i). \qquad (C.2) $$

Once the E-step is performed we can estimate the matrix $\mathbf{\Lambda}$ by

$$ \frac{\partial \chi^2}{\partial \mathbf{\Lambda}} = 0, \qquad (C.3) $$

which gives:

M-step:

$$ \tilde{\mathbf{A}} = \left[ \left( \sum_{i=0}^{N} \mathbf{h}_i \mathbf{h}_i^T \otimes \mathbf{W}_{\mathbf{M}_i} \right) \right]^{-1} \left[ \sum_{i=0}^{N} \mathbf{W}_{\mathbf{M}_i} \mathbf{M}_i \mathbf{h}_i^T \right], \qquad (C.4) $$

where $\tilde{\mathbf{A}}$ is the vector that contains all columns of $\mathbf{\Lambda}$ put end-to-end and the $\otimes$ operator is the Kronecker product. The estimation of the $\mathbf{\Lambda}$ matrix and the orthogonal projections $\mathbf{h}_i$ is done by reiterating the E-step and M-step until $\chi^2$ convergences.

### C.2. Fixing the degeneracies

The free parameters of Eqs. (5) and (22) that we find through orthogonal distance regression contain degenerate degrees of freedom. These must be fixed in order to present an unique an interpretable solution.
### C.2.1. Degeneracies of the extinction fit

The term describing the extinction, $A_{\lambda_{0,i}}\gamma(\lambda)$, is not directly constrained by external observations. As a consequence, Eq. (5) is invariant under several transformations. The first one is

$$A_{\lambda_{0,i}}^{\text{new}} = A_{\lambda_{0,i}} + \sum_{j} h_{i}^{j} c^{j}, \qquad (C.5)$$

$$\alpha_{\lambda}^{j \text{ new}} = \alpha_{\lambda}^{j} - \gamma_{\lambda} c^{j}. \qquad (C.6)$$

We must therefore fix the $c^j$. A natural choice is to impose that the $A_{\lambda_{0,i}}$ are decorrelated with the $h_i^j$. Indeed, any correlation would imply that the extinction term contains information linked to the intrinsic properties described by $h_i^j$. Thus after each E-step, we imposed a correlation of zero between the $A_{\lambda_{0,i}}$ and the $h_i^j$, which amounts to applying the transformations (C.5) and (C.6) taking for the $c^j$

$$\mathbf{c} = [\text{cov}(\tilde{\mathbf{h}})]^{-1} \text{cov}(\tilde{\mathbf{h}}, \mathbf{A}_{\lambda_0}), \qquad (C.7)$$

where $\tilde{\mathbf{h}}$ is the vector that contains the $h_i$, $\text{cov}(\tilde{\mathbf{h}})$ is the covariance matrix of the $\tilde{\mathbf{h}}$, $\text{cov}(\tilde{\mathbf{h}}, \mathbf{A}_{\lambda_0})$ is the vector that contains the covariances between the $h_i^j$ and the $A_{\lambda_{0,i}}$. The second degeneracy to be fixed is the scale of $\gamma_{\lambda}$: only the product $A_{\lambda_{0,i}}\gamma_{\lambda}$ plays a role, and we then impose after each iteration

$$\gamma_{\lambda_0} = 1 \qquad (C.8)$$

by rescaling accordingly $A_{\lambda_{0,i}}$ and $\gamma_{\lambda}$.

Finally, we need a prescription for the mean value of $A_{\lambda_{0,i}}$; we choose it to be centered at zero. This corresponds to the following transformation

$$A_{\lambda_{0,i}}^{\text{new}} = A_{\lambda_{0,i}} - \frac{1}{N} \sum_{i} A_{\lambda_{0,i}}, \qquad (C.9)$$

$$M_{\lambda,0}^{\text{new}} = M_{\lambda,0} + \frac{\gamma_{\lambda}}{N} \sum_{i} A_{\lambda_{0},i}. \qquad (C.10)$$

which leaves the $\chi^2$ invariant and amounts to placing the average spectrum at the average extinction.

### C.2.2. Degeneracies of the global fit

The estimation of the grey parameter, $\Delta M_{\text{grey } i}$, is affected by different degeneracies which can be fixed after the E or M iterations. We choose to fix them after each E-step. The first one arises from the invariance of Eq. (22) under the transformation:

$$\Delta M_{\text{grey } i}^{\text{new}} = \Delta M_{\text{grey } i} + \sum_{j} h_{i}^{j} c^{j}, \qquad (C.11)$$

$$\alpha_{t,\lambda}^{j \text{ new}} = \alpha_{t,\lambda}^{j} - c^{j}. \qquad (C.12)$$

We must therefore fix the $c^j$. Similarly to the prescription made for the extinction fit, a natural choice is to impose that the $\Delta M_{\text{grey } i}$ and the $h_i^j$ are uncorrelated. Indeed, the interpretation is made easier with $\alpha^j$ containing all the information driven by $h_i^j$, including its global impact on magnitudes. We thus impose after

each E-step that the correlation between the $\Delta M_{\text{grey } i}$ and the $h_i^j$ is zero. This is obtained by applying the transformations (C.11) and (C.12) taking for the $c^j$:

$$\mathbf{c} = [\text{cov}(\tilde{\mathbf{h}})]^{-1} \text{cov}(\tilde{\mathbf{h}}, \Delta\mathbf{M}_{\text{grey}}), \qquad (C.13)$$

where $\tilde{\mathbf{h}}_i$ is the vector that contains the $h_i^j$, $\text{cov}(\tilde{\mathbf{h}})$ is the observed covariance matrix computed on the set of $\tilde{\mathbf{h}}_i$ vectors and $\text{cov}(\tilde{\mathbf{h}}, \Delta\mathbf{M}_{\text{grey}})$ is the vector that contains the covariances between $h_i^j$ and $\Delta\mathbf{M}_{\text{grey}}$. Finally, by convention, the $\Delta M_{\text{grey } i}$ are centered on zeros at each step. This is obtained by the following transformation, which leaves the $\chi^2$ invariant:

$$\Delta M_{\text{grey } i}^{\text{new}} = \Delta M_{\text{grey } i} - \frac{1}{N} \sum_{i} \Delta M_{\text{grey } i}, \qquad (C.14)$$

$$M_{t,\lambda,0}^{\text{new}} = M_{t,\lambda,0} + \frac{1}{N} \sum_{i} \Delta M_{\text{grey } i}. \qquad (C.15)$$

## Appendix D: Matrix dispersion estimation

Once the minimum of the $\chi^2$ defined in Eq. (9) is reached, the dispersion matrix is estimated. For this purpose, we use the same method as in Chotard et al. (2011) and which is described in detail in Chotard (2011). The approach is to calculate $\mathbf{D}$ from the observed dispersion of residuals and to subtract the average dispersion due to uncertainties:

$$\mathbf{D} = \frac{1}{N} \sum_{i} (\delta\mathbf{M}_{i} \delta\mathbf{M}_{i}^{T} - \mathbf{C}_{i}) \qquad (D.1)$$

where $\delta\mathbf{M}_{i}$ are the residuals of the model once the minimum of the $\chi^2$ is reached and is defined for a given wavelength as:

$$\delta M_{\lambda,i} = M_{\lambda,i} - M_{\lambda,0} - \sum_{j} q_{i}^{j} \alpha_{\lambda}^{j} - A_{\lambda_{0},i} \gamma_{\lambda} \qquad (D.2)$$

and $\mathbf{C}_i$ is the covariance matrix that accounts for the total propagation of residuals error and is defined as:

$$\mathbf{C}_i = \begin{pmatrix} \cdots & 0 & \cdots \\ \vdots & \sigma_{\lambda i}^2 & \vdots \\ 0 & \cdots & \cdots \end{pmatrix} + (\sigma_{cal}^2 + \sigma_{z}^2) \begin{pmatrix} 1 & \cdots & 1 \\ \vdots & \vdots & \vdots \\ 1 & \cdots & 1 \end{pmatrix} + \alpha \text{cov}(\mathbf{q}_i) \alpha^T \qquad (D.3)$$

where $\alpha$ is the matrix that contains the intrinsic vectors and is defined as:

$$\alpha = (\alpha^1, \alpha^2, \alpha^3, \dots) \qquad (D.4)$$

In order to ensure that the matrix $\mathbf{D}$ is positive definite, the negative eigenvalues of the matrix are set to zero. Once the matrix $\mathbf{D}$ has been calculated, we add it in the expression of the Eq. (9) in order to recalculate the spectral distribution in energy and we iterate the calculations of $\mathbf{D}$ and the computation of the SED, until reaching the maximum of the restricted maximum likelihood (Guy et al. 2010).
### C.2.1. Degeneracies of the extinction fit

The term describing the extinction, $A_{\lambda_{0,i}}\gamma(\lambda)$, is not directly constrained by external observations. As a consequence, Eq. (5) is invariant under several transformations. The first one is

$$A_{\lambda_{0,i}}^{\text{new}} = A_{\lambda_{0,i}} + \sum_{j} h_{i}^{j} c^{j}, \qquad (C.5)$$

$$\alpha_{\lambda}^{j \text{ new}} = \alpha_{\lambda}^{j} - \gamma_{\lambda} c^{j}. \qquad (C.6)$$

We must therefore fix the $c^j$. A natural choice is to impose that the $A_{\lambda_{0,i}}$ are decorrelated with the $h_i^j$. Indeed, any correlation would imply that the extinction term contains information linked to the intrinsic properties described by $h_i^j$. Thus after each E-step, we imposed a correlation of zero between the $A_{\lambda_{0,i}}$ and the $h_i^j$, which amounts to applying the transformations (C.5) and (C.6) taking for the $c^j$

$$\mathbf{c} = [\text{cov}(\tilde{\mathbf{h}})]^{-1} \text{cov}(\tilde{\mathbf{h}}, \mathbf{A}_{\lambda_0}), \qquad (C.7)$$

where $\tilde{\mathbf{h}}$ is the vector that contains the $h_i$, $\text{cov}(\tilde{\mathbf{h}})$ is the covariance matrix of the $\tilde{\mathbf{h}}$, $\text{cov}(\tilde{\mathbf{h}}, \mathbf{A}_{\lambda_0})$ is the vector that contains the covariances between the $h_i^j$ and the $A_{\lambda_{0,i}}$. The second degeneracy to be fixed is the scale of $\gamma_{\lambda}$: only the product $A_{\lambda_{0,i}}\gamma_{\lambda}$ plays a role, and we then impose after each iteration

$$\gamma_{\lambda_0} = 1 \qquad (C.8)$$

by rescaling accordingly $A_{\lambda_{0,i}}$ and $\gamma_{\lambda}$.

Finally, we need a prescription for the mean value of $A_{\lambda_{0,i}}$; we choose it to be centered at zero. This corresponds to the following transformation

$$A_{\lambda_{0,i}}^{\text{new}} = A_{\lambda_{0,i}} - \frac{1}{N} \sum_{i} A_{\lambda_{0,i}}, \qquad (C.9)$$

$$M_{\lambda,0}^{\text{new}} = M_{\lambda,0} + \frac{\gamma_{\lambda}}{N} \sum_{i} A_{\lambda_{0},i}. \qquad (C.10)$$

which leaves the $\chi^2$ invariant and amounts to placing the average spectrum at the average extinction.

### C.2.2. Degeneracies of the global fit

The estimation of the grey parameter, $\Delta M_{\text{grey} i}$, is affected by different degeneracies which can be fixed after the E or M iterations. We choose to fix them after each E-step. The first one arises from the invariance of Eq. (22) under the transformation:

$$\Delta M_{\text{grey} i}^{\text{new}} = \Delta M_{\text{grey} i} + \sum_{j} h_{i}^{j} c^{j}, \qquad (C.11)$$

$$\alpha_{t,\lambda}^{j \text{ new}} = \alpha_{t,\lambda}^{j} - c^{j}. \qquad (C.12)$$

We must therefore fix the $c^j$. Similarly to the prescription made for the extinction fit, a natural choice is to impose that the $\Delta M_{\text{grey} i}$ and the $h_i^j$ are uncorrelated. Indeed, the interpretation is made easier with $\alpha^j$ containing all the information driven by $h_i^j$, including its global impact on magnitudes. We thus impose after

each E-step that the correlation between the $\Delta M_{\text{grey} i}$ and the $h_i^j$ is zero. This is obtained by applying the transformations (C.11) and (C.12) taking for the $c^j$:

$$\mathbf{c} = [\text{cov}(\tilde{\mathbf{h}})]^{-1} \text{cov}(\tilde{\mathbf{h}}, \Delta\mathbf{M}_{\text{grey}}), \qquad (C.13)$$

where $\tilde{\mathbf{h}}_i$ is the vector that contains the $h_i^j$, $\text{cov}(\tilde{\mathbf{h}})$ is the observed covariance matrix computed on the set of $\tilde{\mathbf{h}}_i$ vectors and $\text{cov}(\tilde{\mathbf{h}}, \Delta\mathbf{M}_{\text{grey}})$ is the vector that contains the covariances between $h_i^j$ and $\Delta\mathbf{M}_{\text{grey} i}$. Finally, by convention, the $\Delta M_{\text{grey} i}$ are centered on zeros at each step. This is obtained by the following transformation, which leaves the $\chi^2$ invariant:

$$\Delta M_{\text{grey} i}^{\text{new}} = \Delta M_{\text{grey} i} - \frac{1}{N} \sum_{i} \Delta M_{\text{grey} i}, \qquad (C.14)$$

$$M_{t,\lambda,0}^{\text{new}} = M_{t,\lambda,0} + \frac{1}{N} \sum_{i} \Delta M_{\text{grey} i}. \qquad (C.15)$$

## Appendix D: Matrix dispersion estimation

Once the minimum of the $\chi^2$ defined in Eq. (9) is reached, the dispersion matrix is estimated. For this purpose, we use the same method as in Chotard et al. (2011) and which is described in detail in Chotard (2011). The approach is to calculate **D** from the observed dispersion of residuals and to subtract the average dispersion due to uncertainties:

$$\mathbf{D} = \frac{1}{N} \sum_{i} (\delta\mathbf{M}_{i} \delta\mathbf{M}_{i}^{T} - \mathbf{C}_{i}) \qquad (D.1)$$

where $\delta\mathbf{M}_{i}$ are the residuals of the model once the minimum of the $\chi^2$ is reached and is defined for a given wavelength as:

$$\delta M_{\lambda,i} = M_{\lambda,i} - M_{\lambda,0} - \sum_{j} q_{i}^{j} \alpha_{\lambda}^{j} - A_{\lambda_{0},i} \gamma_{\lambda} \qquad (D.2)$$

and $\mathbf{C}_i$ is the covariance matrix that accounts for the total propagation of residuals error and is defined as:

$$\mathbf{C}_i = \begin{pmatrix} \cdots & 0 & \cdots \\ \vdots & \sigma_{\lambda i}^2 & \vdots \\ 0 & \cdots & \cdots \end{pmatrix} + (\sigma_{cal}^2 + \sigma_{z}^2) \begin{pmatrix} 1 & \cdots & 1 \\ \vdots & \vdots & \vdots \\ 1 & \cdots & 1 \end{pmatrix} + \alpha \text{cov}(\mathbf{q}_i) \alpha^T \qquad (D.3)$$

where $\alpha$ is the matrix that contains the intrinsic vectors and is defined as:

$$\alpha = (\alpha^1, \alpha^2, \alpha^3, \dots) \qquad (D.4)$$

In order to ensure that the matrix **D** is positive definite, the negative eigenvalues of the matrix are set to zero. Once the matrix **D** has been calculated, we add it in the expression of the Eq. (9) in order to recalculate the spectral distribution in energy and we iterate the calculations of **D** and the computation of the SED, until reaching the maximum of the restricted maximum likelihood (Guy et al. 2010).