HYPERION: Hyperspectral Penetrating-type Ellipsoidal Reconstruction for Terahertz Blind Source Separation

Our method, HYPERION, is designed based on Löwner-John ellipsoid (LJE) and the linear mixing model ((Supplementary Note 14)), where the THz data follows a convex geometry structure (see “Methods” for the details). By utilizing the information of the LJE defined as the maximum-volume ellipsoid inscribed in the convex geometry, HYPERION has the efficacy to blindly unmix sources from the data with low data purity [22]. To validate the linear mixing model on the transmitted THz signals, we have measured and characterized mixture tablets composed of D-lactose monohydrate, D-sucrose, and L-tyrosine through a broadband THz time-domain spectroscopy (THz-TDS) system. The three chemicals are selected due to their abundant information within THz spectral bands, including spectral absorption lines and varying absorption spectral trends. The obtained ternary dataset contains the three pure-substance THz spectra (D-lactose monohydrate, D-sucrose, and L-tyrosine) and six mixture spectra composed by each pair of the chemicals with a ratio of 3:7, which is designed to test the data purity criterion. The details of measurement protocol and environment condition are elaborated in Supplementary Note 1.

We compared the THz spectra simulated by the convex combinations of the three chemicals with their corresponding THz measured spectrum with the ternary dataset. All samples are prepared based on the standardized protocol (see “Methods” for the details). As shown in Fig. 1(a), all the simulated spectra are highly aligned with the measured spectra, including the spectral features and the locations of spectral lines. We then project high dimensional data of the ternary dataset into a low dimensional space through principal component analysis (PCA). The PCA evaluation result provides solid evidence that the necessary criteria of holding the linear mixing model are well-preserved (Supplementary Note 2). To further explore the linear mixing model in terms of spectrum similarity and accuracy, we use spectral angle mapper (SAM) and root-mean-square error (RMSE) [23, 24] to evaluate the diminutive differences between simulated spectrum and measured spectrum (Supplementary Note 3). The SAM indicator is based on the projection angle between two spectra in comparison, which explicitly reveals the spectral shape similarity without mixing the influence of spectrum offset. The RMSE indicator is the square root of the quadratic mean of differences between the two spectra. It includes both shape similarity and the effects of the band-wise differences among the measured spectral range. The SAM and RMSE of this ternary dataset are less than 5 degrees and 6 cm^-1 as shown in Fig. 1(b), which demonstrates both the excellent spectral shape similarity as well as the slight deviation between the measured and simulated mixture spectra. Based on the validation of qualitative and quantitative analysis, the assumption of the linear mixing model is experimentally verified. Note that the linear mixing model has been extensively used in remote sensing for reflected (reflecting-type) signals but never for transmitted (penetrating-type) signals as far as we know; thus, the experiment is needed. We then transform the ternary dataset from the THz transmission spectrum to the material absorption spectrum of each measured sample (Supplementary Note 4), which directly presents broadband material information, including resonant frequencies, spectral features, and light-matter interaction levels in the THz range. Additionally, we extract and present material absorption spectrum between 0.2 THz to 1.75 THz since the ASOPS THz-TDS system dynamic range in power spectra of the selected materials approximately decreases to about 5 dB at 1.75 THz that start to mismatch with the linear mixing model. As shown in Fig. 1(c), all unmixed pure substance absorption spectra are highly aligned with the measured ones, including the spectral peak locations and the material absorption profiles. The unmixed D-glucose and L-tyrosine values are slightly deviated from the measured values above 1.5 THz since the signal-to-noise ratio (SNR) of the THz spectroscopy system gradually decreases from 60 dB at 0.3 THz to 5 dB at 1.75 THz.

Nevertheless, at low-SNR spectral regions (1.2 - 1.6 THz), HYPERION still provides great unmixing capability - less than 3.51 cm^-1 deviated absorption coefficient values from the measured ones. For the spectral regimes with a noticeable unmixing difference, for example, the measured values of L-tyrosine absorption peak are 15.35 cm^-1 and 4.78 cm^-1 differences from the unmixed values at 0.949 THz and 1.295 THz, respectively. Due to the overlap of water vapor and material absorption peaks, the linear mixing model is not well-preserved. Since HYPERION takes all the spectral information into account by the LJE instead of specific spectral components, even in nonlinear matter-matter or light-matter interaction spectral regions, the unmixed material absorption spectrum based on HYPERION still demonstrates great spectral accuracy of each source spectral signature but just a relatively higher mismatch in absorption coefficient values.

Despite the nonlinear interaction and the low SNR levels in certain bands, HYPERION offers an accurate unmixing performance – 10.63 degrees with SAM and a 2.54 cm^-1 with RMSE of overall unmixed material absorption spectrum has been achieved. This SAM value presents high spectral shape similarity of unmixed absorption coefficients, preserving rich material spectral profiles and details among a wide frequency range. In another perspective, the extraordinarily low RMSE shows both high unmixing accuracy and spectral similarity of the unmixed absorption coefficient, providing various types of material spectral features for further chemical identification and non-destructive evaluation applications.

HYPERION has several advantages, including support of complex mixture systems, milder data purity requirement, excellent spectral accuracy on unmixed sources, and remarkable noise immunity, due to the simultaneous use of all spectral information by LJE. To further investigate the performance of HYPERION in broad application scope, we design three experimental modules to verify HYPERION characteristics in terms of unmixing capability of complex mixture system, low data purity dataset, and noisy environment condition based on the quinary dataset (see “Methods” for the details). The quinary dataset contains the five pure chemicals (D-glucose, D-lactose monohydrate, L-tyrosine, L-histidine, and D-sucrose) with rich material information in THz spectral range and ten mixture tablets. The five pure chemicals are specifically selected to validate HYPERION performance under different unmixed spectral scenarios - identifying single/multiple spectral peaks and overlaid spectral profiles simultaneously. Within the spectral range of 0.2 - 1.7 THz, L-tyrosine and L-histidine have a single absorption peak at 0.95 and 0.78 THz, respectively. D-lactose monohydrate contains two absorption peaks located at 0.52 THz and 1.3 THz. D-sucrose and D-glucose present the high spectral similarity of their THz absorption spectra, which is used to evaluate ill-conditioned levels of HYPERION. Furthermore, D-lactose is also selected to investigate the unmixing performance of HYPERION on the pharmaceutical application since D-lactose is the common ingredient as an excipient in tablets [25, 26]. The ten mixture tablets are composed of each pair of the five chemicals with a 5:5 ratio, which satisfies the milder requirement of the data purity for a quinary case [27]. Identical experimental setup, tablet preparation, data acquisition, and data preprocessing protocol are well-followed as the previous ternary dataset establishment.

Comparing the unmixed and measured absorption spectra shows that HYPERION has great efficacy on the different spectral scenarios. As shown in Fig. 2, the five unmixed absorption spectra (solid lines) based on the quinary dataset demonstrate a high correlation with the measured absorption spectra (dotted lines) within a broad frequency range. In terms of spectral peak locations, unmixed absorption spectra with varieties of material properties are precisely resolved. The spectral peak deviation of all unmixed absorption spectra is less than 10 GHz, close to the spectral resolution limit of the asynchronized optical sampling (ASOPS) THz-TDS system (see “Methods” for the details) at protocol settings. Additionally, the increasing trends and the characteristics of the absorption spectra are well unmixed since HYPERION utilizes the convex geometry of the whole measured frequency bands for the blind separation instead of specific bands information, such as peak regimes. Those mismatched values between unmixed and measured absorption spectra mainly come from nonlinear light-matter interaction [28], and limited SNR levels at higher frequencies.

To address the ill-conditioned case, HYPERION transforms the data matrix into the preconditioned space, where the similar material absorption spectra can be easily separated (see “Methods” for the details). As shown in Fig. 2, D-glucose and D-sucrose absorption spectra are highly correlated except for the slightly different slope within the 0.2 THz – 1.2 THz frequency range. Under this condition, HYPERION again shows the unmixing capability of a SAM of 9.57 degrees and an RMSE of 2.30 cm^-1, demonstrating its resilience to ill-conditioned scenarios with the complex mixture system. Considering different spectral characteristics among the quinary dataset, HYPERION still reaches 12.15 degrees and 2.88 cm^-1 in SAM and RMSE, respectively. Compared with less complicated material systems (e.g., ternary system), HYPERION maintains a similar accuracy level of unmixed material spectrum based on complex mixture systems (e.g., quinary system). Additionally, HYPERION is featured for its low computation time due to the convexity nature of the LJE (see “Methods” for the details). The computation time of HYPERION on this dataset is less than 3 seconds under a general personal laptop (Intel Core i5 with 8 GB memory).

In most on-site blind source separation scenarios, such as pharmaceutical inspection, biomarker analysis, remote sensing, and chemical identification, the existence and identification of pure substances remain unknown, which generally involve highly mixed chemicals and imply a low data purity. To this end, it is crucial to extract source signals from mixtures with low data purity [22].

Here, we have evaluated the HYPERION based on the same quinary dataset but excluded all pure substances. There is no prerequisite information from all mixture tablets in this testbed, including spectral information, the existence of pure substances, and spectra composition. To address this type of dataset with scarce information, HYPERION fits a simplex to the data convex geometry in the preconditioned space. When the optimal simplex is constructed, the corners of the optimal simplex correspond to the unmixed sources (viewed as vectors) in the preconditioned space (see “Methods” for the details). As shown in Fig. 3, the difference of unmixed absorption spectra among cases with and without pure substances is almost negligible. This is because LJE has been demonstrated to provide comparable performance with the pure substance case [29], as long as the requirement of data purity is satisfied (Supplementary Note 5). The slight difference in performance comes from the measurement noise of the dataset and the number of measured datasets – 15 versus 10 measured tablets. Without the need for pure substance information and its measurement on-site, HYPERION has the potential to expand the THz spectroscopy application scopes, opening up the door to industrial applications in non-invasive sensing, chemical identification, and in vivo biomarker extraction.

To evaluate the efficacy of HYPERION with THz spectroscopy systems in general, we further investigated the noise immunity performance of HYPERION. In this study, we apply additive white Gaussian noise (AWGN) to the measured THz time-domain electric field signals with a 0.001% - 0.1% standard deviation (SD) range (Supplementary Note 10). Moreover, the quinary THz dataset without pure substances is chosen to better evaluate real-world scenarios – low data purity and noisy condition. In those conditions, the data convex geometry will be distorted, resulting in the inaccuracy of unmixed material absorption spectra. To overcome the issue, HYPERION includes a general regularizer in the objective function to accommodate the distortion of convex geometry. With the general regularizer and the regularization parameter, $\lambda$, HYPERION does not force unmixed spectra to form the theoretically required regular structure in the noiseless case. Instead, it encourages the regular structure (see “Methods” section for the details). Based on the design of the regularizer, HYPERION is well-suitable for the penetrating-type THz data, which is normally ill-conditioned and noisy. As shown in Fig. 4(a), the SAM and RMSE of HYPERION unmixing results maintain at low levels while the noise SD is less than 0.1%. In the cases of noisy environment (e.g., top six highest SD conditions), HYPERION performs slightly inferior since the SNR of THz signals in the high-frequency regime drop dramatically at higher noise levels, which causes severe distortion of the data convex geometry and leads to a deteriorating impact of the unmixing performance. To further evaluate the noise immunity capability of HYPERION among unmixed chemicals and frequency ranges, the unmixed absorption spectra of the five pure substances are demonstrated at 0.1% noise SD (Fig. 4(b)-(f)). Under this severe noise condition, the unmixed absorption spectra through HYPERION still show great alignment with the noise-free ground truth in low-frequency bands. Within the 0.2 – 1 THz frequency range, the SAM and RMSE of the D-glucose absorption spectrum are 12.32 degrees and 2.50 cm^-1, respectively. In higher frequency bands, a noticeable difference of every material absorption spectrum starts to show. Although the global trends and features are mostly well-preserved, the high-frequency spectral fluctuations due to the severe additive noise could be an issue for high-precision chemical identification applications.

In sum, HYPERION demonstrates its noise immunity feature with distorted convex geometry datasets. By utilizing this powerful feature, HYPERION not only can fit in a strictly controlled data acquisition environment but also be well suitable for real-world blind detection schemes under a noisy environment.

To further evaluate the unmixing performance of the geometry-based HYPERION, we introduce the commonly used unmixing methods as a comparison: statistic-based nonnegative independent component analysis (nICA), algebra-based nonnegative matrix factorization (NMF), model-based hard modeling factor analysis (HMFA), and geometry-based successive projection algorithm (SPA) [30, 31, 20, 32, 21, 33]. nICA is one of the modified versions of ICA, where it imposes the nonnegative constraint on unmixed sources. With the nonnegative constraint, nICA can unmix the spectra of nonnegative values and lead to better performance than ICA since the fraction of incident radiation absorbed by the material is always nonnegative, leading to the nonnegative values of a material absorption spectrum. NMF is the unmixing method for blind separation problems designed according to the non-convex optimization approaches (Supplementary Note 13). HMFA is the unmixing method based on peak fitting and has demonstrated its efficacy in mid-infrared and THz bands [32, 19]. SPA uses the selection and projection in the vector space to unmix the blind sources. The comparison among HYPERION, nICA, NMF, HMFA, and SPA under different noise SD is shown in Fig. 5. From the unmixed material spectra and the RMSE values, HYPERION has outperformed the THz unmixing methods under a broad noise SD range from 0.001% to 0.1%. As expected, it comes from the three designed properties of HYPERION. First, HYPERION filters out most noise in the affine fitting step, where the data is projected into the lower dimensional convex hull by the singular value decomposition approach. Second, HYPERION uses convex geometry to find the LJE, which can further address the ill-conditioned nature of THz source signals. Third, the embedded general regularizer provides better accommodation to the noisy data by soft fitting for the distortion of convex data geometry in the objective function. As shown in Fig. 5(a, e), although nICA can resolve the spectral absorption peaks, the unmixed material absorption spectra by nICA are severely compressed. The reason is that the nICA method assumes all material absorption spectra are independent in the complex mixture system [31]. However, the summation of the chemical composition (in %) must equal to one, making it impossible for the sources to be statistically independent, and this sum-to-one constraint leads to an assumption of nICA independent assumption. In Fig. 5(b), NMF shows decent unmixing capability because it can converge to the local optimum [34, 35]. However, NMF is relatively sensitive to noisy conditions, probably due to its non-convexity nature. As a result, the unmixing performance of NMF on the THz dataset bounces dramatically at different unmixing trails and noise levels. Consequently, it is challenging to accurately resolve unmixing material spectra under a low SNR spectral regime. As shown in Fig. 5(f), the unmixed absorption spectra through NMF deviate more from the ground truth signal compared to HYPERION while frequency increases. In the case of HMFA (Fig. 5(c, g)), it is not capable of resolving the ill-conditioned data (i.e., D-Glucose and D-Sucrose) nor the spectral trends since the peak fitting approach inherently has difficulties in differentiating overlapped features. As shown in Fig. 5(d, h), SPA has the inferior unmixing performance in the absorption peak regions since the pure-pixel assumption does not hold in those datasets [36]. The convergence time of these methods is also provided in Supplementary Table. 8. Overall, HYPERION delivers superior performance than the four THz unmixing modalities in resolving broadband spectral features, unmixing ill-condition data, reducing noise impact, and preserving spectral details in a noisy environment.

In the sections above, HYPERION has demonstrated superior features on unmixing penetrating-type THz material spectra. Those accurately unmixed spectra can benefit several applications. To demonstrate one of the potential applications, we propose the combination of the raster-scanning THz-TDS system and HYPERION (THz HYPERION) to visualize the “secret recipes” of arbitrary objects at a remote distance. We use the unmixed spectra from the quinary dataset without pure substances and prepare a small test set of 15 chemical tablets with different compositions based on D-lactose monohydrate, D-glucose, and L-tyrosine. The 15 tablets are mounted on the paper board covered by a copper foil, which is utilized to block THz signals encoding the material information of the paper board. Every tablet is measured by the THz-TDS system with the configurations described earlier. The images in visible light and THz band images of the test set are shown in Fig. 6(a) and 6(b), respectively. The unmixed material spectra by the quinary dataset without pure substances are used to restore the chemical compositions of the 15 tablets in the test set by the following convex optimization problem:

	$\displaystyle\min_{\bm{r}_{i}\in\mathbb{R}^{3}}{}$	$\displaystyle~{}\left\|\bm{x}_{i}-\bm{A}\bm{r}_{i}\right\|_{1}$		(1)
	s.t.	$\displaystyle~{}\ \bm{r}_{i}\geq 0,\ \bm{1^{T}}\bm{r}_{i}=1,$

where $i=1,2,\dots,15$. $\bm{x}_{i}$ and $\bm{r}_{i}$ are the measured THz material absorption spectra of the test set and the optimization variables to estimate the corresponding ratio of compositions, respectively. Columns of $\bm{A}$ are the unmixed spectra from the quinary dataset without pure substances. The $L^{1}$ norm is chosen to effectively decrease the influence of the unmixed spectra deviation in the absorption peak regions. Since equation (1) is a convex problem, the solver CVX [37] is adopted for the following demonstration (Supplementary Note 6). The qualitative comparisons between the ground truth and the estimated material composition map are shown in Fig. 6(c, d). As expected, the majority of tablet material compositions are highly aligned with the ground truth except very few tablets (i.e. 10${}^{\text{th}}$, 11${}^{\text{th}}$ and 12${}^{\text{th}}$ tablets). This experimental result indicates that THz HYPERION is applicable for estimating the arbitrary composition of complex mixture systems while the pure material spectra are accurately unmixed. In Fig. 6(e), we further evaluate the unmixed material composition map qualitatively. Most unmixed results present less than 20% inaccuracy of each composition difference compared with ground truth. The inaccuracy of estimating material composition is caused by the number of spectra in the quinary dataset, the high absorption spectrum similarity of D-lactose monohydrate and D-glucose, and the increased system noise in the high-frequency regime. It is worth mentioning that the accuracy of unmixed spectra would be greatly enhanced while the size of the measured dataset is increased. To this extent, the quality of unmixed material signatures can be significantly elevated by extending the measurement area and the pixel number of the multi-material objects under test. Furthermore, different regularizers can be further designed to match with HYPERION for specific application purposes. Overall, our demonstration shows the potential for some blind separation mapping applications, including pharmaceutical analysis, functional imaging, and space exploration.

In complex mixture systems, the material absorption spectra of the overlapped material absorption peaks in the measured THz bands can be quite similar (i.e., D-sucrose and D-glucose). It is the so-called ill-conditioned HU problem, having drawn attention from very recent machine learning literature [42]. As far as we know, the most effective solution for addressing the ill-conditioned HU is based on the Löwner-John ellipsoid (LJE) theory [42], which elegantly exploits the data convex structure (i.e., the convex hull of signature vectors contains the data vectors; cf. equation (2)). However, almost all the HU theories (including LJE) were developed for the reflecting-type signals, probably because the need for the HU techniques mainly comes from the satellite hyperspectral remote sensing, for which the hyperspectral signals are reflected from the objects to the satellite sensors [43]. To apply the LJE theory on the transmitted THz signals, our first task is to reveal the data convex structure of the THz signals.

Let $m_{i}(t),~{}i=1,\dots,n$ be the impulse response function of the $i^{\text{th}}$ material that the input THz signal $x(t)$ penetrates through, and let $\otimes$ denote the convolution operator. The transmitted THz signal measured at time $t$ can then be modeled as

$$y(t)=x(t)\otimes m_{1}(t)\otimes\cdots\otimes m_{n}(t).$$

By taking the Fourier transform on both sides, we then have

$$Y(f)=X(f)~{}\!M_{1}(f)\cdots M_{n}(f),$$

where $(Y,X,M_{i})$ are the Fourier transforms of $(y,x,m_{i})$, respectively, and $f$ is the frequency index. If there are $k$ frequency samples $f_{1},\dots,f_{k}$ in the THz spectral regions, then $(Y,X,M_{i})$ can be considered as $k$-dimensional vectors (e.g., $Y=[Y(f_{1}),\dots,Y(f_{k})]^{T}$). Let us define the standardized data as $Z(f)=\log\left(\frac{|Y(f)|}{|X(f)|}\right)=\sum_{i=1}^{n}{\log(|M_{i}(f)|)}$. Also, the materials $M_{1},\dots,M_{n}$ may not be distinct, we assume that there are $q$ distinct materials in the set $\{M_{1},\dots,M_{n}\}$, and let $N(i)$ be the set of indices corresponding to the $i^{\text{th}}$ distinct material in the set $\{M_{1},\dots,M_{n}\}$. By defining $H_{i}(f)=\prod_{j\in N(i)}M_{j}(f)$, we can further simplify the expression as $Z(f)=\sum_{i=1}^{q}{\log(|H_{i}(f)|)}$. Here, we neglect the energy loss between the interfaces of materials since it is relatively small compared to energy loss in the lossy medium (Supplementary Note 4). According to the derived approximation of material absorption coefficients (Supplementary Note 4), we have $\log(|H_{i}(f)|)=\frac{1}{2}\alpha_{i}(f)d_{i},~{}i=1,\dots,q$, thereby leading to the standardized data representation

$$Z(f)=\sum_{i=1}^{q}\frac{1}{2}\alpha_{i}(f)d_{i},$$

where $\alpha_{i}(f)$ is the material absorption coefficient of the $i^{\text{th}}$ distinct material at frequency $f$, and $d_{i}$ is the penetration depth of the $i^{\text{th}}$ distinct material.

Naturally, we define the hyperspectral signature of the $i^{\text{th}}$ material as

$$\bm{s}_{i}\triangleq[\alpha_{i}(f_{1}),\dots,\alpha_{i}(f_{k})]^{T}\in\mathbb{R}^{k},~{}i=1,\dots,q.$$

Let $d_{i}^{n}$ be the penetration depth of the $i^{\text{th}}$ distinct material in the $n^{\text{th}}$ sample data. It is elegant to observe that if we have the information of the value of $l_{n}\triangleq\frac{1}{2}(d_{1}^{n}+\dots+d_{q}^{n})$, simply normalizing the $n^{\text{th}}$ standardized data $Z_{n}\triangleq[~{}\!Z_{n}(f_{1}),\dots,Z_{n}(f_{k})~{}\!]^{T}$ by $l_{n}$ can reveal the desired convex structure, i.e., $\frac{Z_{n}}{l_{n}}\in\textrm{conv}\{{\bm{s}}_{1},\dots,{\bm{s}}_{q}\}$. If the incident angle is zero, $l_{n}$ is nothing but half of the thickness of the $n^{\text{th}}$ sample. Therefore, in the proposed HYPERION algorithm, the $\ell$ samples $Z_{1},\dots,Z_{\ell}$ are allowed to have different thickness $2l_{1},\dots,2l_{\ell}$. Based on different thickness of samples, the normalized/standardized data ${\bm{x}}_{n}=\left[\frac{Z_{n}(f_{1})}{l_{n}},\dots,\frac{Z_{n}(f_{k})}{l_{n}}\right]^{T}$ have the desired data convex structure

$${\bm{x}}_{1},\dots,{\bm{x}}_{\ell}\in\textrm{conv}\{{\bm{s}}_{1},\dots,{\bm{s}}_{q}\}.$$

(2)

The aim is to recover the THz signatures ${\bm{s}}_{1},\dots,{\bm{s}}_{q}$ from the preprocessed THz transmitted signals ${\bm{x}}_{1},\dots,{\bm{x}}_{\ell}$ for material identification.

Due to the convex structure in equation (2), HYPERION can use the information of LJE, defined as the maximum-volume ellipsoid inscribed in the standardized/normalized penetrating-type THz signals as shown in Fig. 7, for the source separation with mild requirement of data purity (Supplementary Note 5).

To be precise, the standardized/normalized transmitted THz data matrix ${\bm{X}}$ can be explicitly written as

$${\bm{X}}\triangleq[{\bm{x}}_{1},\dots,{\bm{x}}_{\ell}]={\bm{P}}{\bm{T}}\bm{\Sigma}\in\mathbb{R}^{k\times\ell},$$

where ${\bm{P}}\in\mathbb{R}^{k\times q}$ with $[{\bm{P}}]_{ij}\triangleq\alpha_{j}(f_{i})$, ${\bm{T}}\in\mathbb{R}^{q\times\ell}$ with $[{\bm{T}}]_{ij}\triangleq\frac{1}{2}d_{i}^{j}$, and $\bm{\Sigma}\in\mathbb{R}^{\ell\times\ell}$ is the diagonal matrix with the $i^{\text{th}}$ diagonal entry being $\frac{1}{l_{i}}$. Then, computing the LJE of the THz signals ${\bm{X}}$ can be proven to be a convex optimization problem [44] as below:

	$\displaystyle(\bm{F}^{\star},\bm{c}^{\star})=$	$\displaystyle\arg\max_{{\bm{F}}\in\mathbb{S}^{q-1}_{++},~{}\!{\bm{c}}\in\mathbb{R}^{q-1}}~{}\log\det({\bm{F}})$		(3)
		$\displaystyle{\rm s.t.}~{}~{}~{}~{}~{}~{}\lVert\bm{F}\bm{b}_{i}\rVert\leq h_{i}-\bm{b}_{i}^{T}\bm{c},~{}\forall i=1,\dots,H,$

where $\mathbb{S}^{p-1}_{++}$ is the positive semidefinite (PSD) cone, the halfspace parameters $\{({\bm{b}}_{1},h_{1}),\dots,({\bm{b}}_{H},h_{H})\}$ come from the $\cal H$-polytope representation of the transmitted THz data ${\bm{X}}$, i.e., $\mathcal{H}({\bm{X}})=\mathcal{H}({\bm{P}}{\bm{T}}\bm{\Sigma})\equiv\mathcal{H}(\textrm{conv}\{\bm{x}_{1},\dots,\bm{x}_{\ell}\})=\bigcap_{i=1}^{H}~{}\!\{\bm{x}\mid\bm{b}_{i}^{T}\bm{x}\leq h_{i}\}$ [45], and the optimal argument $(\bm{F}^{\star},\bm{c}^{\star})$ gives the desired maximum-volume ellipsoid inscribed in the THz data ${\bm{X}}\equiv\mathcal{H}({\bm{X}})$ (i.e., LJE) as

$$\mathcal{E}(\bm{F}^{\star},\bm{c}^{\star})\triangleq\{\bm{F}^{\star}\bm{\alpha}+\bm{c}^{\star}\mid\|\bm{\alpha}\|\leq 1\}\subseteq\mathbb{R}^{q-1}.$$

Note that we have assumed w.l.o.g. that the dimension of $\textrm{conv}\{\bm{x}_{1},\dots,\bm{x}_{\ell}\}$ is no greater than $q-1$, as implied by equation (2). The solving details for equation (3) are discussed in the Supplementary Note 7.

As aforementioned, THz data is ill-conditioned. As reported in [42], the information embedded in the LJE $\mathcal{E}({\bm{F}}^{\star},{\bm{c}}^{\star})$ is critical in preconditioning the data ${\bm{X}}$; specifically, LJE yields the preconditioned data

$$\mathcal{C}({\bm{X}})=({\bm{F}}^{\star})^{\dagger}\left({\bm{P}}{\bm{T}}\bm{\Sigma}-{\bm{c}}^{\star}\bm{1}_{\ell}^{T}\right)\in\mathbb{R}^{(q-1)\times\ell},$$

whose corresponding (preconditioned) THz signatures are easily verified to be $\widetilde{{\bm{S}}}=({\bm{F}}^{\star})^{\dagger}\left([{\bm{s}}_{1},\dots,{\bm{s}}_{q}]-{\bm{c}}^{\star}\bm{1}_{q}^{T}\right)$; note that the precondition operator $\mathcal{C}(\cdot)$ is applied columnwisely as $\mathcal{C}({\bm{v}})=({\bm{F}}^{\star})^{\dagger}({\bm{v}}-{\bm{c}}^{\star})$ for any given (column) vector ${\bm{v}}$. Remarkably, estimating the preconditioned THz signatures $\widetilde{{\bm{S}}}$ is much more friendly because, according to the LJE theory [42, Theorem 1], columns of $\widetilde{{\bm{S}}}$ will form a regular simplex centered at the origin whenever the data purity $\gamma>\frac{1}{\sqrt{q-1}}$. As $\widetilde{{\bm{S}}}$ is obtained, the THz signatures can be simply recovered as $[{\bm{s}}_{1},\dots,{\bm{s}}_{q}]=({\bm{F}}^{\star})^{\dagger}\widetilde{{\bm{S}}}+{\bm{c}}^{\star}\bm{1}_{q}^{T}$. Therefore, we can focus on estimating $\widetilde{{\bm{S}}}$ from the preconditioned THz data $\mathcal{C}({\bm{X}})$ next.

According to the regular simplex structure, $\widetilde{{\bm{S}}}$ can be characterized as [42, Equation (11)]

$$\widetilde{{\bm{S}}}=\alpha{\bm{U}}^{T}{\bm{S}}_{0},$$

(4)

where ${\bm{S}}_{0}$ forms any unit-volume regular simplex centered at the origin $\bm{0}_{q-1}$, with closed-form expression available in [42, Proposition 2]; the unitary matrix ${\bm{U}}$ (i.e., ${\bm{U}}^{T}{\bm{U}}=\bm{I}_{q-1}$) and the scalar $\alpha>0$ (with closed-form expression available in [42, Equation (12)]) are for rotating and scaling ${\bm{S}}_{0}$ to fit $\widetilde{{\bm{S}}}$.

As the transmitted THz data is not just ill-conditioned but also noisy, the regular structure may not be strictly satisfied. Therefore, unlike [42], we did not force the regular structure. Instead, we just use such information to design a regularizer $\phi(\widetilde{{\bm{S}}})\triangleq\left\|\widetilde{{\bm{S}}}-\alpha\bm{U}^{T}{\bm{S}}_{0}\right\|_{F}^{2}$, by minimizing which we just encourage the regular structure (rather than forcing it).

To finish the design of the HU criterion for the transmitted THz data, we need to design the data-fitting term. To this end, we need the following lemma, whose proof is given in Supplementary Note 8.

By Lemma 1, we naturally design the data fitting term as $\left\|\mathcal{C}({\bm{X}})-\widetilde{{\bm{S}}}{\bm{T}}\bm{\Sigma}\right\|_{F}^{2}$, which, together with the regularizer $\phi$, yields the following HU criterion for transmitted THz signals:

	$\displaystyle\min_{\widetilde{{\bm{S}}},~{}\!{\bm{T}},\bm{\Sigma},~{}\!\bm{U}}~{}$	$\displaystyle~{}\left\|\mathcal{C}({\bm{X}})-\widetilde{{\bm{S}}}{\bm{T}}\bm{\Sigma}\right\|_{F}^{2}+\lambda\phi(\widetilde{{\bm{S}}})+I_{\textrm{unitary}}({\bm{U}})$		(5)
	s.t.	$\displaystyle~{}{\bm{T}}\geq\bm{0},~{}\bm{\Sigma}\geq\bm{0},~{}\bm{1}_{q}^{T}{\bm{T}}\bm{\Sigma}=\bm{1}_{\ell}^{T},$

where $\lambda:=1$ is the regularization parameter, and $I_{\textrm{unitary}}(\cdot)$ is the indicator function of the set of unitary matrices. An algorithm for solving equation (5) is provided in Supplementary Note 9. Once equation (5) is solved, the column vectors of the optimal solution $\widetilde{{\bm{S}}}^{\star}$ then serves as the estimates of the THz signatures in the preconditioned space as shown in Fig. 7. Based on the problem formulation and the solving detail, the algorithm, termed HYperspectral Penetrating-type Ellipsoidal ReconstructION (HYPERION), has been completed. Remarkably, HYPERION does not use any information about the pattern of resonant peaks of the signatures and is designed under a fully unsupervised setting.

In the asynchronized optical sampling (ASOPS) THz-TDS system (Menlo TERA ASOPS, Menlo Systems, Germany) as shown in Fig. 8, two asynchronized Er-doped fiber femtosecond lasers are fed into an InGaAs/InAlAs THz photoconductive antenna emitter and an LT-InGaAs/InAlAs THz detector. The power, bandwidth, and dynamic range of the THz-TDS system are up to 60 µW, less than 4.5 THz, and greater than 80 dB. The generated THz radiation from the THz photoconductive antenna emitter then consecutively travels through two convex THz lenses with a focal length of 10 cm. The first convex lens is used for THz beam collimation, and the second convex lens is to focus the THz radiations on the tablets. The diameter of the focused beam is approximately 1.5 mm of full width at half maximum (FWHM). After penetrating through tablets, the transmitted THz waves contained material absorption information in THz range, then traveling through two identical convex THz lenses to the THz photoconductive antenna detector. THz photoconductive antenna detector is used to retrieve the time-resolved THz electric field and convert the electric field signal to the photocurrent. The connected transimpedance amplifier (TIA) then amplifies the induced photocurrent to voltage signals. The bandwidth and the gain of the TIA are 1.8 MHz and 10⁶ (V/A), respectively. Each measurement contains a 100 ps time-domain trace with a 5 fs temporal resolution. To further increase the signal-to-noise ratio (SNR) of the system, we average 1,000 time-domain traces, which can effectively decrease the time-domain additive white Gaussian noise to 0.015% according to the law of large numbers (Supplementary Note 10). With the setting above, the system offers a dynamic range of over 65 dB from 0.1 THz to 3 THz.

In this study, all tablets were composed of D-glucose, D-lactose monohydrate, L-histidine, L-tyrosine, and D-sucrose powders. D-glucose, D-lactose monohydrate, and D-sucrose powders were from Merck & Co., Inc. (Kenilworth, NJ, USA). L-histidine, L-tyrosine powders were from Acros Organics (Geel, Belgium). Each chemical was ground and mixed uniformly by the ball mill, ensuring that the Mie scattering is greatly reduced in the measured spectral characteristics. Followed by the pulverization, the chemical powders were poured into a tablet die made by a high-speed steel (HSS) for subsequent high-pressure compression. Powders were compressed by a hydraulic press (Specac, Orpington, U.K.) to form 3 mm-thick tablets (total mass: 0.3 g; thickness: $3.05\pm 0.02$ mm; see Supplementary Table 6 for the details) under the pressure of 1,000 kg/cm² for 15 seconds. Tablets were then mounted on a 0.6 mm thick polylactic acid (PLA) plate. The copper foil is covered around the tablets to block THz signals encoding the PLA plate material information. Additionally, the copper foil prevents the low-frequency diffraction effect caused by the edge of the PLA plate.

Each spectrum is obtained by fast Fourier transform (FFT) of a 100 picoseconds (ps) THz time-domain trace based on an asynchronous optical sampling (ASOPS) THz-TDS system (see Supplementary Note 11). To increase the system dynamic range, 1,000 spectra of each mixture are measured and averaged, providing more than 60 dB dynamic range at 0.3 THz and 25 dB dynamic range at 1.75 THz. All the measurements are conducted at room temperature (23°C), one atmosphere (atm) pressure, and 60% humidity. Under this condition, the dataset contains water vapor absorption lines at 0.56 THz, 0.75 THz, 0.99 THz, 1.10 THz, 1.16 THz, 1.21 THz, 1.23 THz, 1.41 THz, 1.60 THz, 1.66 THz, and 1.72 THz within 0.2 - 1.75 THz frequency range (Supplementary Fig. 9), which contribute nearly 2.3% of the nonlinear region in the measured frequency bands.

Whether the dataset follows the linear mixing model is essential since HYPERION utilizes the convex geometry of the model to unmix the material signatures. To validate the linear mixing model upon the dataset, principle component analysis (PCA) is one of the most suitable tools since data information can still be reserved after projecting to a lower-dimensional space [46]. To be precise, we project the non-pure mixture signatures to the optimal 2 dimensional (2D) space, which is constructed by the two eigenvectors corresponding to the two largest eigenvalues. Additionally, only non-pure mixture points are selected for the calculation of the 2D space since the assumption of pure substances is not held. After obtaining the optimal 2D plane, the data points of pure substances are further projected to the optimal 2D plane. In the evaluation of the linear mixing model, the quinary dataset with pure substances is selected due to the material variety. As shown in Supplementary Fig. 10, the high-dimensional data points are projected onto the two-dimensional plane. Except for the extreme cases, most of the data points in every subplot form the triangle in the projected 2D space, where three points lie at the vertex and three points fall at the center of the side. The misaligned data points are caused by the inaccurate estimation of the two-dimensional plane, which is caused by the measurement noise and nonlinear light-matter interaction. Based on the projected data points, the measured data points are highly aligned with the linear mixing model. Additionally, the triangles involving D-glucose and D-lactose are narrower since the signatures of the two chemicals are the ill-conditioned case. It is worthwhile to mention that the alignment accuracy between the model and measured data may be further enhanced by modeling the nonlinear matter-matter interaction into the linear mixing model.

Spectral angle mapper (SAM) is typically used to evaluate the shape similarity between two spectra. The definition is as a following equation:

$$\theta(s_{1},s_{2})=\cos^{-1}\left(\frac{s_{1}\cdot s_{2}}{\left\|s_{1}\right\|_{2}\cdot\left\|s_{2}\right\|_{2}}\right),$$

(6)

where $s_{1}$ and $s_{2}$ are spectra.
Root mean square error (RMSE) does not only focus on the shape similarity but the band-wise differences. RMSE is defined as a following equation:

$$d(s_{1},s_{2})=\sqrt{\frac{\left\|s_{1}-s_{2}\right\|_{2}^{2}}{n}},$$

(7)

where $n$ is the dimension of the spectrum.

The THz-TDS system provides a time-domain waveform that profiles the interactions between THz radiation and measured tablets. The reference waveform $S_{\text{ref}}(t)$ and the tablet waveform $S_{\text{tablet}}(t)$, where $t$ is the time index, refer to the measured waveform without and with a tablet, respectively. After obtaining reference waveform $S_{\text{ref}}(t)$ and tablet waveform $S_{\text{tablet}}(t)$, the fast Fourier transform (FFT) was applied to convert time-domain waveforms to the frequency-domain reference spectrum $S_{\text{ref}}(\omega)$ and tablet spectrum $S_{\text{tablet}}(\omega)$, where $\omega$ is the angular frequency. To further extract the material absorption spectrum ($\alpha(\omega)$) of the tablet, the thickness of the tablet is an important parameter for the derivation. The thickness of the tablet is measured by the electric ruler, which provides a 0.1 mm resolution. Compared to the thickness of the tablet (i.e.,  3 mm), this resolution of 0.1 mm introduces approximately 3% inaccuracy in the measurement. The material absorption spectra of the tablets are derived from the Fresnel equations and the THz wave propagation model based on homogeneous and planar materials. Furthermore, according to this derivation, the thickness measurement error of 3% will only introduce approximately 3% error in material absorption spectra. Since the particle size of the tablet powder is vastly smaller than the wavelength of the THz waves (Supplementary Note 16) Additionally, the thickness of the tablet is significantly larger than the effective wavelength of the THz signal in the 0.2 - 1.75 THz frequency range such that the Fabry-Pérot effect [17] does not introduce drastic damping to the measured spectrum.

Under this condition, we do not include the scattering influence and the Fabry-Pérot effect in this model. The Fresnel equations describe the transmission and the reflection of a THz wave traveling through an interface, which are based on the complex refractive index of the material, $\tilde{n}(\omega)=n(\omega)-j\kappa(\omega)$, where $n(\omega)$ is the real refractive index and $\kappa(\omega)$ represents the extinction coefficient which is proportional to the material absorption coefficient, $\kappa(\omega)=\left[\alpha(\omega)\cdot\lambda)\right]/4\pi$, where $\lambda$ is the wavelength. The material absorption coefficients describe the power loss when the THz waves travel through the material. Considering a typical incident THz wave propagates a material with thickness $d$ (Supplementary Fig. 11), the Fresnel equation [17] at an interface can be represented as

$$t_{ab}(\omega)=\frac{2\tilde{n}_{a}(\omega)}{\tilde{n}_{a}(\omega)+\tilde{n}_{b}(\omega)},$$

$$r_{ab}(\omega)=\frac{\tilde{n}_{a}(\omega)-\tilde{n}_{b}(\omega)}{\tilde{n}_{a}(\omega)+\tilde{n}_{b}(\omega)},$$

where $t_{ab}(\omega)$ is the transmission coefficient of a THz wave from the region $a$ (which is air in our experiment) to region $b$ (which is the tablet in our experiment) and $r_{ab}(\omega)$ is the reflection coefficient of a THz wave at the $a$-$b$ interface. Additionally, $\tilde{n}_{a}(\omega)$ and $\tilde{n}_{b}(\omega)$) represent the complex refractive index of region $a$ and $b$, respectively. The THz wave propagation constant, $p_{b}(\omega,d)$, is governed by

$$p_{b}(\omega,d)=\exp\left[\frac{-j\tilde{n}_{b}(\omega)\omega d}{c}\right],$$

where $c$ is the speed of light. Consequently, the equation for the received THz tablet signal in the frequency domain can be represented as

	$\displaystyle S_{\text{tablet}}(\omega)$	$\displaystyle=S_{\text{ref}}(\omega)\cdot t_{ab}\cdot t_{ba}\cdot p_{b}(\omega,d)\quad\textbf{(The Fresnel loss is included in $t_{ab}$ and $t_{ba}$.)}$
		$\displaystyle=S_{\text{ref}}(\omega)\cdot t_{ab}\cdot t_{ba}\cdot\exp\left[\frac{-j\tilde{n}_{b}(\omega)\omega d}{c}\right]$
		$\displaystyle=S_{\text{ref}}(\omega)\cdot t_{ab}\cdot t_{ba}\cdot\exp\left\{\frac{-j[n_{b}(\omega)-j\kappa_{b}(\omega)]\omega d}{c}\right\}$
		$\displaystyle=S_{\text{ref}}(\omega)\cdot t_{ab}\cdot t_{ba}\cdot\exp\left[\frac{\kappa_{b}(\omega)\omega d}{c}\right]\cdot\exp\left[\frac{-jn_{b}(\omega)\omega d}{c}\right].$

The reference signal normalizes the tablet signal, and its amplitude can be obtained from equation (6)

$$\left|\frac{S_{\text{tablet}}(\omega)}{S_{\text{ref}}(\omega)}\right|=\left|t_{ab}\cdot t_{ba}\cdot\exp\left[\frac{\kappa_{b}(\omega)\omega d}{c}\right]\right|\triangleq\left|H(\omega)\right|.$$

(8)

The extinction coefficient can be derived by taking the natural logarithm operation of equation (8)

$$\ln(\left|H(\omega)\right|)=\frac{\kappa_{b}(\omega)\omega d}{c}+\ln(\left|t_{ab}\cdot t_{ba}\right|),$$

$$\kappa_{b}(\omega)=\frac{c\left[\ln(\left|H(\omega)\right|)-\ln(\left|t_{ab}\cdot t_{ba}\right|)\right]}{\omega d}.$$

Therefore, the material absorption coefficient can be represented as

$$\alpha^{{}^{\prime}}(\omega)=\frac{4\pi\kappa(\omega)}{\lambda}=\frac{2\cdot\left[\ln(\left|H(\omega)\right|)-\ln(\left|t_{ab}\cdot t_{ba}\right|)\right]}{d}.$$

(9)

Based on the derivation, to acquire the material absorption spectrum, we have first to calculate the transmission coefficients, $t_{ab}$ and $t_{ba}$. However, the interface power loss is relatively small since the thickness of the sample is greatly larger than the THz wavelength, and the tablet surface is smooth and flat compared to the THz wavelength. This relatively small interface power loss results in the high transmission coefficients of $t_{ab}$ and $t_{ba}$. Thus, we can neglected the second term of the right-hand side in equation (9) to form the equation (10) since the multiplication of $t_{ab}$ and $t_{ba}$ approaches to $1$.

$$\alpha^{{}^{\prime}}(\omega)=\frac{2\cdot\ln(\left|H(\omega)\right|)}{d}=\frac{2\cdot\ln\left(\left|\frac{S_{\text{tablet}}(\omega)}{S_{\text{ref}}(\omega)}\right|\right)}{d}.$$

(10)

To further match the notation in blind source separation fields, the equation (10) can be expressed in ordinary frequency, $f$, by the equation of $\omega=2\pi f$:

$$\alpha^{{}^{\prime}}(\omega)=\alpha^{{}^{\prime}}(2\pi f)\triangleq\alpha(f).$$

(11)

The aforementioned problem (i.e., recovering ${\bm{s}}_{1},\dots,{\bm{s}}_{q}$ from their mixtures ${\bm{x}}_{1},\dots,{\bm{x}}_{\ell}$) has been studied in the machine learning literature; for example, in the successive projection algorithm (SPA) [47], the problem is solved under one additional assumption called separability assumption, i.e.,

$${\bm{s}}_{1},\dots,{\bm{s}}_{q}\in\{{\bm{x}}_{1},\dots,{\bm{x}}_{\ell}\},$$

meaning that, for each material, there is a pure observation solely composed of such material. Under the separability assumption, one can show that $\textrm{conv}\{{\bm{x}}_{1},\dots,{\bm{x}}_{\ell}\}=\textrm{conv}\{{\bm{s}}_{1},\dots,{\bm{s}}_{q}\}$, and hence identifying the signatures is nothing but seeking the vertices of the observable geometry $\textrm{conv}\{{\bm{x}}_{1},\dots,{\bm{x}}_{\ell}\}$, which is polynomial-time solvable. Remarkably, by following the definition of data purity $\gamma\in(0,1]$, defined in [27], the separability assumption holds if, and only if, $\gamma=1$ (the highest data purity). So, the separability assumption is quite restricted, and such assumption is often violated in practical hyperspectral applications [27].

In fact, the separability assumption may not hold in our THz application, as a tablet typically involves highly mixed chemicals, implying a low data purity. In recent mathematics literature, the LJE theory has been proven to work well even with very low data purity; specifically, the signatures can be perfectly recovered under a much milder assumption $\gamma>\frac{1}{\sqrt{q-1}}$ [22]. Also, as aforementioned, the LJE criterion is theoretically and experimentally proven to be robust against the ill-conditioning of the signatures, hence quite suitable for our THz applications. Therefore, in this paper, the convex LJE criterion is applied to the transmitted (penetrating-type) THz signals for the first time, to be detailed later.

We remark that there is other non-convex criterion that can also perform HU under the mild requirement of $\gamma>\frac{1}{\sqrt{q-1}}$ [48, 27], but such criterion cannot achieve the theoretical bound $\frac{1}{\sqrt{q-1}}$ in practice due to its non-convexity [22]. By contrast, the LJE criterion has the same theoretical bound, and in the meanwhile, it is experimentally shown to be able to achieve the bound due to its convexity nature [22].

To restore the chemical compositions, the equation (12) should be solved for every tablet.

	$\displaystyle\min_{\bm{r}_{i}\in\mathbb{R}^{3}}{}$	$\displaystyle~{}\left\|\bm{x}_{i}-\bm{A}\bm{r}_{i}\right\|_{1}$		(12)
	s.t.	$\displaystyle~{}\ \bm{r}_{i}\succeq 0,\ \bm{1^{T}}\bm{r}_{i}=1.$

Since the equation (12) is convex, it can be solved by the solver CVX [37] as the following pseudo-code:

	$\displaystyle{\rm cvx\underline{~{}~{}}begin}$
	$\displaystyle~{}~{}~{}~{}~{}~{}{\rm variable~{}r(3)~{}nonnegative;}$
	$\displaystyle~{}~{}~{}~{}~{}~{}{\rm minimize(~{}norm(x-A*r,1)~{});}$
	$\displaystyle~{}~{}~{}~{}~{}~{}{\rm subject~{}to}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm sum(r)==1;}$
	$\displaystyle{\rm cvx\underline{~{}~{}}end}$

One of the core in HYPERION is the LJE problem:

	$\displaystyle(\bm{F}^{\star},\bm{c}^{\star})=$	$\displaystyle\arg\max_{{\bm{F}}\in\mathbb{S}^{q-1}_{++},~{}\!{\bm{c}}\in\mathbb{R}^{q-1}}~{}\log\det({\bm{F}})$		(13)
		$\displaystyle{\rm s.t.}~{}~{}~{}~{}~{}~{}\lVert\bm{F}\bm{b}_{i}\rVert\leq h_{i}-\bm{b}_{i}^{T}\bm{c},~{}\forall i=1,\dots,H,$

where $\mathbb{S}^{p-1}_{++}$ is the positive semidefinite (PSD) cone, the halfspace parameters $\{({\bm{b}}_{1},h_{1}),\dots,({\bm{b}}_{H},h_{H})\}$ come from the $\cal H$-polytope representation of the data ${\bm{X}}$ and the optimal argument $(\bm{F}^{\star},\bm{c}^{\star})$ gives the desired maximum-volume ellipsoid inscribed in the penetrating-type THz data.

We remark that there is an algorithm [42] that can be adapted to solve equation (13) exactly, but it is mainly designed to handle the case with million-scale constraints $H$ induced by NASA’s benchmark hyperspectral data [49]. In typical THz applications, such as measuring the chemicals in a given tablet, we just have a few samples, and hence we do not need to use such a sophisticated method. Instead, the general-purpose convex software CVX [37] is effective and fast when the problem size is not too large, and hence is suitable to be adopted to solve equation (13).

To adopt CVX, we need to implement the constraints and the objective function of equation (13). The constraints of equation (13) are the well-known Lorentz cone constraints “$\lVert\bm{F}\bm{b}_{i}\rVert\leq h_{i}-\bm{b}_{i}^{T}\bm{c}$”, theoretically equivalent to require the LJE to lie within the penetrating-type data convex hull $\mathcal{H}({\bm{P}}{\bm{T}}\bm{\Sigma})\equiv\mathcal{H}(\textrm{conv}\{\bm{x}_{1},\dots,\bm{x}_{\ell}\})$. As for the objective function, for PSD matrix ${\bm{F}}$, we observe that

$$\log\det({\bm{F}})\propto\det({\bm{F}})\propto(\det({\bm{F}}))^{\frac{1}{q-1}}.$$

So, the objective function can be implemented by the command “${\rm det\underline{~{}~{}}rootn(\cdot)}$” supported by the CVX package. Solving the LJE problem by the general-purpose convex software CVX normally reuqires some prior knowledge to transform the problem into the supporting functions. To this end, we provide the pseudo-code for the better understanding. The pseudo-code using CVX is as following:

	$\displaystyle{\rm cvx\underline{~{}~{}}begin}$
	$\displaystyle~{}~{}~{}~{}~{}~{}{\rm variable~{}F(q-1,q-1)~{}symmetric;}$
	$\displaystyle~{}~{}~{}~{}~{}~{}{\rm variable~{}c(q-1);}$
	$\displaystyle~{}~{}~{}~{}~{}~{}{\rm maximize(~{}det\underline{~{}~{}}rootn(F)~{});}$
	$\displaystyle~{}~{}~{}~{}~{}~{}{\rm subject~{}to}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm for~{}i=1:H}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm norm(~{}F*b(i,:),~{}2~{})+b(i,:)^{\prime}*c<=h(i);}$
	$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm end}$
	$\displaystyle{\rm cvx\underline{~{}~{}}end}$

It shows that CVX can handle a non-standard form of convex programming.

Following the definition of precondition operator $\mathcal{C}$, we have

	$\displaystyle\mathcal{C}({\bm{X}})$	$\displaystyle\triangleq({\bm{F}}^{\star})^{\dagger}\left({\bm{P}}{\bm{T}}\bm{\Sigma}-{\bm{c}}^{\star}\bm{1}_{\ell}^{T}\right)$
		$\displaystyle=({\bm{F}}^{\star})^{\dagger}{\bm{P}}{\bm{T}}\bm{\Sigma}-({\bm{F}}^{\star})^{\dagger}{\bm{c}}^{\star}\bm{1}_{\ell}^{T}$
		$\displaystyle=({\bm{F}}^{\star})^{\dagger}{\bm{P}}{\bm{T}}\bm{\Sigma}-({\bm{F}}^{\star})^{\dagger}{\bm{c}}^{\star}\bm{1}_{q}^{T}{\bm{T}}\bm{\Sigma}$
		$\displaystyle=({\bm{F}}^{\star})^{\dagger}[{\bm{s}}_{1},\dots,{\bm{s}}_{q}]{\bm{T}}\bm{\Sigma}-({\bm{F}}^{\star})^{\dagger}{\bm{c}}^{\star}\bm{1}_{q}^{T}{\bm{T}}\bm{\Sigma}$
		$\displaystyle=\left[({\bm{F}}^{\star})^{\dagger}\left([{\bm{s}}_{1},\dots,{\bm{s}}_{q}]-{\bm{c}}^{\star}\bm{1}_{q}^{T}\right)\right]{\bm{T}}\bm{\Sigma}$
		$\displaystyle=\widetilde{{\bm{S}}}{\bm{T}}\bm{\Sigma},$

where we have used the observation that each column of the matrix ${\bm{T}}\bm{\Sigma}$ actually provides a set of convex combination coefficients [44], implying the equation “$\bm{1}_{q}^{T}{\bm{T}}\bm{\Sigma}=\bm{1}_{\ell}^{T}$”. $\blacksquare$

As there are four variables in equation (14),

	$\displaystyle\min_{\widetilde{{\bm{S}}},~{}\!{\bm{T}},\bm{\Sigma},~{}\!\bm{U}}{}$	$\displaystyle~{}\left\|\mathcal{C}({\bm{X}})-\widetilde{{\bm{S}}}{\bm{T}}\bm{\Sigma}\right\|_{F}^{2}+\lambda\phi(\widetilde{{\bm{S}}})+I_{\textrm{unitary}}({\bm{U}})$		(14)
	s.t.	$\displaystyle~{}{\bm{T}}\geq\bm{0},~{}\bm{\Sigma}\geq\bm{0},~{}\bm{1}_{q}^{T}{\bm{T}}\bm{\Sigma}=\bm{1}_{\ell}^{T},$

solving it directly using the alternating optimization may not be effective. Considering that what we are really interested in is the THz signatures $\widetilde{{\bm{S}}}$, we do not need to separately optimize all the variables, allowing us to introduce the trick of changes of variables to reduce the number of block variables. Specifically, by letting $\widetilde{{\bm{T}}}={\bm{T}}\bm{\Sigma}$, equation (14) can be reformulated as

	$\displaystyle\min_{\widetilde{{\bm{S}}},~{}\!\widetilde{{\bm{T}}},~{}\!\bm{U}}{}$	$\displaystyle~{}\left\|\mathcal{C}({\bm{X}})-\widetilde{{\bm{S}}}\widetilde{{\bm{T}}}\right\|_{F}^{2}+\lambda\phi(\widetilde{{\bm{S}}})$		(15)
	s.t.	$\displaystyle~{}{\widetilde{{\bm{T}}}\geq\bm{0},~{}\bm{1}_{q}^{T}\widetilde{{\bm{T}}}=\bm{1}_{\ell}^{T},}~{}\bm{U}^{T}\bm{U}=\bm{I}_{q-1}.$

After reducing the number of block variables, it is natural to solve equation (15) using alternating optimization that alternatively updates $\widetilde{{\bm{S}}},~{}\!\widetilde{{\bm{T}}},~{}\!\bm{U}$, as detailed next.

The subproblem for solving $\widetilde{{\bm{S}}}$ (with $\widetilde{{\bm{T}}},~{}\!\bm{U}$ fixed) can be explicitly written as

$\displaystyle\widetilde{{\bm{S}}}^{\star}=\arg\min_{\widetilde{{\bm{S}}}}\left\|\mathcal{C}({\bm{X}})-\widetilde{{\bm{S}}}\widetilde{{\bm{T}}}\right\|_{F}^{2}+\lambda\left\|\widetilde{{\bm{S}}}-\alpha{\bm{U}}^{T}{\bm{S}}_{0}\right\|_{F}^{2}.$

(16)

By defining $\widetilde{\bm{m}}\triangleq\textrm{vec}(\widetilde{{\bm{S}}})$, ${\bm{v}}\triangleq[\textrm{vec}(\mathcal{C}({\bm{X}}))^{T},~{}\sqrt{\lambda}(\textrm{vec}(\alpha{\bm{U}}^{T}{\bm{S}}_{0}))^{T}]^{T}$ and $\widetilde{{\bm{P}}}\triangleq[\widetilde{{\bm{T}}}\otimes\bm{I}_{q-1},~{}\sqrt{\lambda}\bm{I}_{q(q-1)}]^{T}$, equation (16) can be simplified as $\textrm{vec}(\widetilde{{\bm{S}}}^{\star})\triangleq\widetilde{\bm{m}}^{\star}=\arg\min_{\widetilde{\bm{m}}}~{}\|\widetilde{{\bm{P}}}\widetilde{\bm{m}}-{\bm{v}}\|_{2}^{2}$, which closed-form solution can be derived as

	$\displaystyle\widetilde{\bm{m}}^{\star}$	$\displaystyle=(\widetilde{{\bm{P}}}^{T}\widetilde{{\bm{P}}})^{-1}(\widetilde{{\bm{P}}}^{T}{\bm{v}})$
		$\displaystyle=\{[\widetilde{{\bm{T}}}\widetilde{{\bm{T}}}^{T}+\lambda\bm{I}_{q}]^{-1}\otimes\bm{I}_{q-1}\}\textrm{vec}(\mathcal{C}({\bm{X}})\widetilde{{\bm{T}}}^{T}+\lambda\alpha{\bm{U}}^{T}{\bm{S}}_{0})$
		$\displaystyle=\textrm{vec}\left(\bm{I}_{q-1}~{}\![\mathcal{C}({\bm{X}})\widetilde{{\bm{T}}}^{T}+\lambda\alpha{\bm{U}}^{T}{\bm{S}}_{0}]~{}\![\widetilde{{\bm{T}}}\widetilde{{\bm{T}}}^{T}+\lambda\bm{I}_{q}]^{-T}\right),$

where $\otimes$ and $\textrm{vec}(\cdot)$ denote Kronecker product and vectorization operator, respectively. Therefore, we have $\widetilde{{\bm{S}}}^{\star}=[\mathcal{C}({\bm{X}})\widetilde{{\bm{T}}}^{T}+\lambda\alpha{\bm{U}}^{T}{\bm{S}}_{0}][\widetilde{{\bm{T}}}\widetilde{{\bm{T}}}^{T}+\lambda\bm{I}_{q}]^{-T}$.

Next, note that the subproblem for solving $\widetilde{{\bm{T}}}$ (with $\widetilde{{\bm{S}}},~{}\!\bm{U}$ fixed) is

$\displaystyle\widetilde{{\bm{T}}}^{\star}=\arg\min_{\widetilde{{\bm{T}}}\geq\bm{0},~{}\bm{1}_{q}^{T}\widetilde{{\bm{T}}}=\bm{1}_{\ell}^{T}}{}$

$\displaystyle~{}\left\|\mathcal{C}({\bm{X}})-\widetilde{{\bm{S}}}\widetilde{{\bm{T}}}\right\|_{F}^{2},$

which can be easily solved as the well-known fully-constrained least-squares problem [50, 51]. Finally, the subproblem for solving ${\bm{U}}$ (with $\widetilde{{\bm{S}}},~{}\!\widetilde{{\bm{T}}}$ fixed) is a non-convex optimization problem

$\displaystyle{\bm{U}}^{\star}=\arg\min_{\bm{U}^{T}\bm{U}=\bm{I}_{q-1}}\left\|\widetilde{{\bm{S}}}-\alpha{\bm{U}}^{T}{\bm{S}}_{0}\right\|_{F}^{2},$

(17)

but can be elegantly solved using the following inequalities for square matrices ${\bm{A}},{\bm{B}}$, i.e., $\textrm{trace}({\bm{A}})\leq\|{\bm{A}}\|_{*}=\|{\bm{A}}\|_{S^{1}}$ and $\|{\bm{A}}{\bm{B}}\|_{S^{1}}\leq\|{\bm{A}}\|_{S^{p}}\|{\bm{B}}\|_{S^{q}}$, where $\frac{1}{p}+\frac{1}{q}=1$ ($1\leq p,q\leq\infty$), $\|\cdot\|_{*}$ denotes nuclear norm, and $\|\cdot\|_{S^{p}}$ denotes Schatten $p$-norm. Since ${\bm{U}}$ is unitary, the objective function of equation (17) can be simplified as $\textrm{trace}({\bm{U}}\widetilde{{\bm{S}}}{\bm{S}}_{0}^{T})$, whose upper bound is derivable using the inequalities as $\textrm{trace}({\bm{U}}\widetilde{{\bm{S}}}{\bm{S}}_{0}^{T})\leq\|{\bm{U}}\widetilde{{\bm{S}}}{\bm{S}}_{0}^{T}\|_{S^{1}}\leq\|{\bm{U}}\|_{S^{\infty}}\|\widetilde{{\bm{S}}}{\bm{S}}_{0}^{T}\|_{S^{1}}=\|\widetilde{{\bm{S}}}{\bm{S}}_{0}^{T}\|_{S^{1}}$, and is achievable by ${\bm{U}}^{\star}={\bm{D}}_{1}{\bm{D}}_{2}^{T}$, where ${\bm{D}}_{1},{\bm{D}}_{2}$ are obtained from the singular value decomposition of $\widetilde{{\bm{S}}}{\bm{S}}_{0}^{T}={\bm{D}}_{2}\bm{\Sigma}^{\prime}{\bm{D}}_{1}^{T}$.

The algorithm, termed HYperspectral Penetrating-type Ellipsoidal ReconstructION (HYPERION), has been completed. Remarkably, HYPERION does not use any information about the pattern of resonant peaks and is designed under a fully blind setting.

The noise in our THz-TDS system corresponds to additive white Gaussian noise (AWGN), whose standard deviation in the single time-domain trace is approximately 0.48% of the peak-to-peak amplitude as shown in Supplementary Fig. 12. According to the law of large numbers, the standard deviation of AWGN is linearly proportional to $\sqrt{n}$, where $n$ is the number of time traces. In this regard, we applied the AWGN to our dataset to simulate the different levels of noise conditions. We also present this noise amplitude in terms of SNR by adopting the equation (18) from the speech signal processing field.

$\displaystyle\text{SNR}=10\log_{10}\left(\frac{P(S)-P(N)}{P(N)}\right),$

(18)

where $P(S)$ and $P(N)$ are the power of signal and noise, respectively.

The conventional THz time-domain spectroscopy (THz-TDS) system is essential for many THz applications, such as remote sensing, material characterization, imaging, and defect inspection. As THz-TDS systems are capable of sampling ultrafast THz signals in the time domain, THz-TDS systems can provide more information, such as time delay, spectral phase, and spectral amplitude [52]. A conventional THz-TDS system utilizes a mechanical delay-line stage to introduce a time delay between two split femtosecond laser beams: the so-called pump beam and probe beam [53]. The pump beam is coupled to the THz photoconductive antenna emitter to generate the THz radiations. The probe beam is fed to the THz photoconductive antenna detector to gate received THz signals [54, 52]. With the accumulated time delay between two beams, a THz-TDS system can profile the THz time-domain signal, where the effective sampling frequency is typically tens of THz, which is the reciprocal of the time delay [53]. However, the sampling time and the quality of the time-domain signal are severely limited by the mechanical delay line stage condition, such as maximum moving velocity, position resolution, and operational stability. To address the limitation, the asynchronized optical sampling (ASOPS) THz-TDS system is developed to increase the sampling frequency. Additionally, the ASOPS THz-TDS system is more stable due to the elimination of the mechanical delay line stage interference [55]. By replacing the mechanical delay line stage with the two asynchronized femtosecond lasers, the time delay between the probe beam and the pump beam can be accumulated in a short time. More specifically, the sampling time of a single time trace, $t_{\text{sample}}$, can be decreased down to within millisecond timescale. Additionally, $t_{\text{sample}}$ is related to the repetition rate of the lasers as

$\displaystyle t_{\text{sample}}=\frac{\frac{1}{f^{A}_{\text{rep}}}\times\frac{1}{f^{A}_{\text{rep}}}}{\frac{1}{f^{A}_{\text{rep}}}-\frac{1}{f^{A}_{\text{rep}}+\Delta f}}\cong\frac{1}{\Delta f},$

where $f^{A}_{\text{rep}}$ and $\Delta f$ are the repetition rate of the pump beam laser and repetition rate difference between the two femtosecond lasers. In our ASOPS THz-TDS system, the repetition rate of the pump beam laser is 100 MHz. To achieve the high time-domain sampling resolution and the fast sampling time simultaneously, the repetition rate difference is set as 10 Hz, delivering 5 fs sampling resolution and 100 ms sampling time. However, under this configuration, the size of the full single-time trace is considerably large, requiring high data transfer rate and large storage space. To prevent those issues, we have only taken 100 ps segment of the single time trace and set the sampling rate of the DAQ card as 20 MHz. Accordingly, our ASOPS THz-TDS system provides 0.01 THz resolution in the spectrum and 5 fs time-domain sampling resolution.

The detectable bands of the spectra are between 0.2 THz and 1.75 THz due to physical and hardware limitations [56, 57]; the THz spectral signal below 0.2 THz suffers from the physical and instrumental limitations, such as photoconductive material carrier lifetime, bandwidth, and sensing area of the THz photoconductive antennas; the THz spectral signal above 1.75 THz decays to the system noise level due to the absorption of the tablets. Consequently, we have only extracted THz bands between 0.2 THz and 1.75 THz as our valid spectra for the further unmixing experiments.

The nonnegative matrix factorization (NMF) is one of the commonly used methods to recover material signatures in the blind sense. To address the non-convexity of NMF, the result of NMF is taken from one of the 10 trials under different initialization based on the best root-mean-squared error (RMSE). The RMSE, $D$, is defined as:

$\displaystyle D=\frac{\left\|\bm{X}-\bm{WH}\right\|_{F}}{\sqrt{nm}},$

where $\mathbf{X}\in\mathbb{R}^{n\times m}$, $\mathbf{W}\in\mathbb{R}^{n\times k}$, and $\mathbf{H}\in\mathbb{R}^{k\times m}$ are the target, signatures matrix, and abundance matrix, respectively. In the experiment, $k$ equals to 3 and 5 in the ternary and quinary cases, respectively. Additionally, we adopt the naive NMF in equation (19) since validity of the commonly-used regularizers, such as $L^{2}$ norm or $L^{1}$ norm, have not yet been verified with reasonable physical meaning. The alternating least square method is used to solve this NMF optimization problem:

	$\displaystyle\min_{{\bm{W}},~{}\!{\bm{H}}}{}$	$\displaystyle~{}\left\|{\bm{X}}-{\bm{W}}{\bm{H}}\right\|_{F}$		(19)
	s.t.	$\displaystyle~{}{{\bm{W}}\geq\bm{0},~{}{\bm{H}}\geq\bm{0}.}$

Since the sum-to-one criterion derived in the “Methods” section is based on the assumption where interface power loss is relatively small compared to material absorption. To validate the assumption in practical cases, we have measured the spectra of three mixture tablets with different mixing approaches as shown in Supplementary Fig. 14. The three tablets are composed of the same amount of substances. As shown in Supplementary Fig. 15, the deviation among three measured spectra by different mixing approaches is only less than 2 dB. Based on the experiment, it concludes that the spectrum is not severely affected by the mixing approaches.

In addition to the mixing approaches, the dynamic range of the THz-TDS system affects the accuracy of linear mixing model on the material absorption spectra since the largest detectable absorption coefficients $\alpha_{\text{max}}(f)$ are determined by the THz-TDS system dynamic range [58]. More specifically, the performance of HYPERION will be decreased when the material absorption coefficients are larger than the $\alpha_{\text{max}}(f)$. To this sense, it is important to compare the $\alpha_{\text{max}}(f)$ of our ASOPS THz-TDS system to the pure substance absorption spectra. As shown in Supplementary Fig. 13, in the range between 1.3 THz and 1.5 THz, the material absorption coefficients of D-lactose monohydrate and D-glucose excess the $\alpha_{\text{max}}(f)$. Except this region, the material absorption spectra are less than the $\alpha_{\text{max}}(f)$. Although the measured absorption coefficients are not accurate in this region, HYPERION can still deliver the acceptable unmixing performance due to use of the material absorption information from all measured frequency bands.

We compared the material composition estimation by the material absorption spectra unmixed from different methods, including HYPERION, NMF, HMFA, nICA, and SPA. The configuration of the test set is same as in the “Application” section. As shown in Supplementary Table. 1, HYPERION delivers superior performance of estimating material composition since HYPERION unmixs the accurate material absorption spectra. HMFA and nICA have the much inferior performance on material estimation since the unmixed material absorption spectra are much suppressed in magnitude. The unmixed material spectra by the quinary dataset without pure substances are used to restore the chemical compositions of the 15 tablets in the test set by the convex optimization problem in equation (20).

	$\displaystyle\min_{\bm{r}_{i}\in\mathbb{R}^{3}}~{}$	$\displaystyle~{}\left\|\bm{x}_{i}-\bm{A}\bm{r}_{i}\right\|_{1}$		(20)
	s.t.	$\displaystyle~{}\ \bm{r}_{i}\geq 0,\ \bm{1^{T}}\bm{r}_{i}=1,$		(21)

Here, we will investigate the effect of the scattering on the unmixed spectra results. Since the effect of scattering is proportional to the grain size of the tablet [59], an optical microscopy with 1.25 numerical aperture (NA) is utilized to measure the grain sizes of different tablets. In the measurement, all chemicals follow the same grounding approach (see “Method” for the details) and are placed on the standard microscopy slide. For each chemical, the grain size is determined by taking the average of the grain size of 7 randomly selected grains. As shown in Supplementary Table. 2, the grain sizes of all chemicals fall in 1 µm range. Thus, the scattering effect in different chemicals falls in the same scale. To this sense, we can model the measured material absorption as the addition of true material absorption and the scattering effect as in equation (22).

$\displaystyle u_{msr}(f)=u_{abs}(f)+u_{sct}(f),$

(22)

where $u_{msr}$, $u_{abs}$, and $u_{sct}$ are the measured material absorption, material absorption and the scattering effect, respectively. Since the grain sizes of different chemicals fall in the same range, $u_{sct}(f)$ is considerably identical for different chemicals. In HYPERION, this constant shift $u_{sct}(f)$ will be calibrated in the affine fitting step and added to the unmixed spectrum in the transformation between the preconditioned space and original space. Thus, the scattering effect will not affect the unmixing performance of HYPERION when the grain sizes of different chemical fall in the same scale.