\section{Computational methods}

The main computational task required to implement our approach is the solution of the \(\ell_1\) optimization problem in eq.~\ref{eq:bpdn}. From the many algorithms available for this purpose, we rely on the spectral projected gradient \(\ell_1\) (SPGL1) algorithm developed by van~den~Berg and Friedlander~\cite{Berg2008} and their freely-available implementation.  

For all compressed sensing calculations in this paper, the change-of-basis matrix between the sparse basis and the measurement basis is given by the DCT matrix whose elements are given by
%
\begin{align}
\label{eq:dct}
P_{ij} = \sqrt{\dfrac{2}{N}} \cos \left[ \frac{\pi}{N} \left(i-1\right) \left(j - \frac{1}{2}\right) \right] \ ,
\end{align}
%
with the first row multiplied by an extra factor of \(1/\sqrt{2}\) to guarantee orthogonality.

For the numerical calculations we avoid explicitly constructing the Kronecker product of \(P\) with itself and instead perform all matrix multiplications in the SPGL1 algorithm directly in terms of \(P\).  This latter approach has much smaller memory requirements and numerical costs, ensuring that the compressed sensing process itself is rapid and not a bottleneck in our procedure.  The condition \(PAP^T = B\) is satisfied up to a relative error of \(10^{-7}\) in the Frobenius norm (vectorial 2-norm).

In order to perform the undersampling required for our compressed sensing calculations, first the complete Hessians were calculated, then they were converted to the measurement basis, and finally they were randomly sampled by column.  Quantum mechanical Hessians were obtained with the QChem 4.2~\cite{Shao_2014} software package, using density functional theory with the B3LYP exchange-correlation functional~\cite{Becke_1993} and the 6-31G* basis set.  Molecular mechanics Hessians were calculated using the the MM3 force field~\cite{Allinger_1989} and the open-source package Tinker 6.2.