Non-negative matrix factorization for 2D-XAS images of lithium ion batteries

The data acquired in Nakamura et al [9] were provided by the coauthors of this paper. They two-dimensionally observed the change in the Co-K-edge x-ray absorption spectra of Li_x CoO₂ in the model composite electrode. The XAS measurements at the Co K edge were carried out at the BL01B1 and BL28XU beamlines of SPring-8, Japan. Figures 1 (A) and (B) show two-dimensional maps of the PTE obtained from the data used in this study. An area of about 1100 × 3060 μm was observed with a spatial resolution of approximately 2 μm in an energy range from 7700 to 7740 eV at energy steps of 0.4 eV. The model electrode was charged to 4.2 V at current rates of 9 mAcm⁻² (figure 1 (B)) and 12 mAcm⁻² (figure 1 (A)). The procedure for fabricating the model electrode and the details of the experiments are given in the papers [9, 29]. The two 2D-XAS datasets we analyzed consisted of 100 frames (corresponding to an energy range from 7700 to 7740 eV at energy steps of 0.4 eV) of 300 × 300 pixels (corresponding to area of about 600 × 600 μm), which were regions straddling the boundary between the electrochemically active and inactive regions, as shown by the red squares in figures 1(A) and (B).

In the preprocessing, we normalized each Co K-edge x-ray absorption spectrum at each pixel point to be zero at 7700 eV and unity at 7740 eV. Then, we reshaped the 2D array of pixels into a 1D column of pixels for each frame and thereby obtained an m by n spectrum matrix X ₀. In this study, m = 90000 and n = 100. As described in the Introduction, x-ray absorption spectra acquired from the model electrode are very similar and vary only a little from one position to another. To enhance the differences in the absorption spectra, we subtracted a certain reference standard spectrum from the spectrum matrix X ₀:

$$\begin{eqnarray}&&{\boldsymbol{X}}={{\boldsymbol{X}}}_{0}-{\bf{1}}{\boldsymbol{b}},\end{eqnarray} \tag{ 1 }$$

where 1 is an m-dimensional column vector in which all the elements are 1, and b is an n-dimensional row vector storing the reference standard spectrum. The matrix X obtained by this subtraction is called the difference spectrum matrix. The spectrum matrix X ₀ is non-negative, whereas the difference spectrum matrix X takes signed values due to the subtraction in equation (1). Here, we used the x-ray absorption spectrum of Li_0.5CoO₂ reference material (figure 1(C)) as the reference standard spectrum to detect sensitively spectral changes due to charging and discharging, since Li_0.5CoO₂ corresponds to the fully charged state [9]. The reference spectrum in figure 1(C) was separately obtained from the powder sample of Li_0.5CoO₂ through the transmission spectrum measurement with a high intensity beam over a sufficient measurement time.

We assume that the difference spectrum matrix X can be factorized into a real nonnegative matrix ${\boldsymbol{W}}\in {{\mathbb{R}}}_{+}^{m\times r}$ and a real signed matrix ${\boldsymbol{H}}\in {{\mathbb{R}}}^{n\times r}$ whose rank is equal to $r\lt \min (m,n)$. Under this assumption, we introduce the following generative model for X ,

$$\begin{eqnarray}&&{\boldsymbol{X}}={\boldsymbol{W}}{{\boldsymbol{H}}}^{T}+{\boldsymbol{\varepsilon }},\end{eqnarray} \tag{ 2 }$$

where all elements of ${\boldsymbol{\varepsilon }}\in {{\mathbb{R}}}^{m\times n}$ are independent and identically distributed (i.i.d.) Gaussian noise. The k-th column vector of W , denoted as ${{\boldsymbol{w}}}_{k}\in {{\mathbb{R}}}_{+}^{m\times 1}$, represents the k-th spatial domain (k ∈ {1, ⋯ ,r}), and the k-th column vector of H , denoted as ${{\boldsymbol{h}}}_{k}\in {{\mathbb{R}}}^{n\times 1}$, represents the k-th difference spectrum corresponding to the k-th spatial domain (k ∈ {1, ⋯ ,r}). To represent the difference spectrum taking signed values due to the subtraction in equation (1), each element of h _k is defined as a signed real value. Using these column and row vectors, the generative model of equation (2) can be rewritten as

$$\begin{eqnarray}&&{\boldsymbol{X}}=\sum _{k=1}^{r}{{\boldsymbol{w}}}_{k}{{\boldsymbol{h}}}_{k}^{T}+{\boldsymbol{\varepsilon }}.\end{eqnarray} \tag{ 3 }$$

Under the i.i.d. Gaussian noise assumption in equation (2), W and H can be estimated by minimizing the following objective function,

$$\begin{eqnarray}&&J=\parallel {\boldsymbol{X}}-{\boldsymbol{W}}{{\boldsymbol{H}}}^{T}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 4 }$$

where ∥X∥₂ denotes the Frobenius norm.

We seek to minimize the objective function (4) subject to column orthogonality on W by using the hierarchical alternating least squares (HALS) algorithm for the ONMF proposed by Kimura et al [24]. We summarize the HALS algorithm for one-sided ONMF below.

The HALS algorithm minimizes the objective function defined in equation (4) alternately with respect to each of the column vectors w _k and h _k in W and H while keeping the other vectors fixed [30]. Each alternating minimization with respect to w _k and h _k can be equivalently formulated as a minimization of the following objective function,

$$\begin{eqnarray}&&{J}_{k}=\parallel {{\boldsymbol{X}}}^{(k)}-{{\boldsymbol{w}}}_{k}{{\boldsymbol{h}}}_{k}^{T}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\boldsymbol{X}}}^{(k)}={\boldsymbol{X}}-\sum _{{k}^{{\prime} }(\ne k)}{{\boldsymbol{w}}}_{{k}^{{\prime} }}{{\boldsymbol{h}}}_{{k}^{{\prime} }}^{T},\end{eqnarray} \tag{ 6 }$$

where X ^(k) denotes the residual between the data X and the factorization result excluding the k-th factor ${{\boldsymbol{w}}}_{k}{{\boldsymbol{h}}}_{k}^{T}$.

Kimura et al formulated an optimization problem that minimizes J_k subject to the following constraint for imposing column orthogonality on one side's nonnegative factor matrix [24]:

$$\begin{eqnarray}&&{{\boldsymbol{W}}}^{(k)T}{{\boldsymbol{w}}}_{k}=\sum _{{k}^{{\prime} }(\ne k)}{{\boldsymbol{w}}}_{{k}^{{\prime} }}^{T}{{\boldsymbol{w}}}_{k}=0,\end{eqnarray} \tag{ 7 }$$

where W ^(k) denotes the sum of all column vectors in W excluding the k-th column vector w _k . To solve this optimization problem, they introduced the following Lagrange function with Lagrange multiplier λ_k ,

$$\begin{eqnarray}&&{L}_{k}=\parallel {{\boldsymbol{X}}}^{(k)}-{{\boldsymbol{w}}}_{k}{{\boldsymbol{h}}}_{k}^{T}{\parallel }_{2}^{2}+{\lambda }_{k}({{\boldsymbol{W}}}^{(k)T}{{\boldsymbol{w}}}_{k}).\end{eqnarray} \tag{ 8 }$$

Then, they found the stationary points at which the first partial derivatives of L_k with respect to w _k and h _k are zero and obtained an update rule of the HALS algorithm for ONMF. The value of the Lagrange multiplier λ_k can be uniquely determined under the column-orthogonality assumption for W [31].

The one-sided ONMF we use here imposes column orthogonality on the nonnegative W , but relaxes the non-negativity restriction on H . To represent the difference spectrum taking signed values, each element of h _k is defined as a real signed value. This is a point of difference from the HALS algorithm proposed by Kimura et al In our case, the HALS algorithm for the one-sided ONMF is given by the following alternating update equation:

$$\begin{eqnarray}&&{{\boldsymbol{w}}}_{k}\leftarrow {\left[{{\boldsymbol{X}}}^{(k)}{{\boldsymbol{h}}}_{k}-\displaystyle \frac{{{\boldsymbol{W}}}^{(k)T}{{\boldsymbol{X}}}^{(k)}{{\boldsymbol{h}}}_{k}}{{{\boldsymbol{W}}}^{(k)T}{{\boldsymbol{W}}}^{(k)}}{{\boldsymbol{W}}}^{(k)}\right]}_{+},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\boldsymbol{h}}}_{k}\leftarrow {{{\boldsymbol{X}}}^{(k)}}^{T}{{\boldsymbol{w}}}_{k},\end{eqnarray} \tag{ 10 }$$

where [ · ]₊ denotes half-wave rectification ${[x]}_{+}=\max (x,0)$ resulting in a nonnegative w _k . Note that there is no such half-wave rectification on h _k , which is different from the HALS algorithm proposed by Kimura et al [24]. λ_k can be set to $\tfrac{{{\boldsymbol{W}}}^{(k)T}{{\boldsymbol{X}}}^{(k)}{{\boldsymbol{h}}}_{k}}{{{\boldsymbol{W}}}^{(k)T}{{\boldsymbol{W}}}^{(k)}}$ as in Kimura et al [24]. Algorithm 1 is the HALS algorithm for the one-sided ONMF.

In the preprocessing, we obtain difference spectra by subtracting the reference standard spectrum from 2D-XAS data. After the preprocessing, we apply the one-sided ONMF to the difference spectra matrix X as described above. Thus, H represents background-subtracted spectra of the extracted domains. To obtain original x-ray absorption spectra of the extracted domains, we calculate W ^T X ₀. Because the one-sided ONMF imposes column orthogonality on one side's nonnegative matrix W representing individual domains, W ^T W ≃ I is satisfied, and thus calculating W ^T X ₀ can approximately give the x-ray absorption spectra of the individual domains. This calculation corresponds to a weighted average of the x-ray absorption spectra for each domain.

Algorithm 1. HALS algorithm for one-sided ONMF

1: Initialize ${\boldsymbol{W}}$ and ${\boldsymbol{H}}$ randomly.

2: ${\boldsymbol{U}}={\boldsymbol{W}}{{\bf{1}}}_{r}$

3: white Convergence criterion is not satisfied do

4: ${\boldsymbol{A}}={\boldsymbol{X}}{\boldsymbol{H}}$

5: ${\boldsymbol{B}}={{\boldsymbol{H}}}^{T}{\boldsymbol{H}}$

6: for k = 1 to r do

7: ${{\boldsymbol{W}}}^{(k)}={\boldsymbol{U}}-{{\boldsymbol{w}}}_{k}$

8: $u={{\boldsymbol{A}}}_{:k}-{\boldsymbol{W}}{{\boldsymbol{B}}}_{:k}+{{\boldsymbol{B}}}_{{kk}}{{\boldsymbol{w}}}_{k}$

9: ${{\boldsymbol{w}}}_{k}={\left[{\boldsymbol{u}}-\tfrac{{{\boldsymbol{W}}}^{(k)T}{\boldsymbol{u}}}{{{\boldsymbol{W}}}^{(k)T}{{\boldsymbol{W}}}^{(k)}}{{\boldsymbol{W}}}^{(k)}\right]}_{+}$

10: ${{\boldsymbol{w}}}_{k}={{\boldsymbol{w}}}_{k}/\parallel {{\boldsymbol{w}}}_{k}{\parallel }^{2}$

11: ${\boldsymbol{U}}={{\boldsymbol{W}}}^{(k)}+{{\boldsymbol{w}}}_{k}$

12: end for

13: ${\boldsymbol{C}}={{\boldsymbol{W}}}^{T}{\boldsymbol{X}}$

14: ${\boldsymbol{D}}={{\boldsymbol{W}}}^{T}{\boldsymbol{W}}$

15: for k = 1 to r do

16: ${{\boldsymbol{h}}}_{k}={{\boldsymbol{C}}}_{:k}-{\boldsymbol{H}}{{\boldsymbol{D}}}_{:k}+{{\boldsymbol{D}}}_{{kk}}{{\boldsymbol{h}}}_{k}$

17: end for

18: end while

We seek to estimate the number of spatial domains with different spectral profiles by applying model selection to the one-sided ONMF. Here, we can use bi-cross-validation [25] to determine the rank of the two factor matrices from the data X. Bi-cross-validation is applicable to a wide class of matrix factorization methods [25].

The data ${\boldsymbol{X}}\in {{\mathbb{R}}}^{m\times n}$ are separated into four blocks consisting of X _A , X _B , X _C and X _D , as follows:

$$\begin{eqnarray}{\boldsymbol{X}}=\left[\begin{array}{cc}{{\boldsymbol{X}}}_{A} & {{\boldsymbol{X}}}_{B}\\ {{\boldsymbol{X}}}_{C} & {{\boldsymbol{X}}}_{D}\end{array}\right].\end{eqnarray} \tag{ 11 }$$

Let X _D be factorized into W _D and H _D by applying the HALS algorithm described above. The estimated value of the denoised X _D becomes ${\hat{{\boldsymbol{X}}}}_{D}={{\boldsymbol{W}}}_{D}{{\boldsymbol{H}}}_{D}^{T}$. As reported in [25], the generalization error can be measured from a residual defined by ${\hat{{\boldsymbol{\varepsilon }}}}_{A}={{\boldsymbol{X}}}_{A}-{{\boldsymbol{X}}}_{B}{\hat{{{\boldsymbol{X}}}_{D}}}^{+}{{\boldsymbol{X}}}_{C}$, where ${\hat{{{\boldsymbol{X}}}_{D}}}^{+}$ is the pseudo-inverse of $\hat{{{\boldsymbol{X}}}_{D}}$.

The detailed procedure for calculating the cross-validation error is as follows. Every 20th row and column vectors are selected from X . The partial blocks X _A , X _B , X _C and X _D are organized with these picked vectors and the remaining row and column vectors. Then, X _D is factorized into W _D and H _D by the HALS described above, and the generalization error $\parallel {\hat{{\boldsymbol{\varepsilon }}}}_{A}{\parallel }_{2}^{2}$ is calculated. This process is repeated 20 times until all the row and column vectors have been picked, and the average of $\parallel {\hat{{\boldsymbol{\varepsilon }}}}_{A}{\parallel }_{2}^{2}$ is obtained. Hereafter, the average of $\parallel {\hat{{\boldsymbol{\varepsilon }}}}_{A}{\parallel }_{2}^{2}$ is called the bi-cross-validation error. Then, the rank of the two factor matrices was determined by searching for the minimum point of the bi-cross-validation error.

We generated synthetic data as ground truth on which we could objectively assess the performance of our method.

We manually created eight 400 by 400 binary arrays representing eight separate domains (figure 2(A)) and stored them as a 160 000 by 8 binary matrix W _g by reshaping the 2D binary array of each of these domains into a 1D binary array. We synthesized eight different spectra on a 100 by 8 matrix H _g , in which each consisted of a randomly shifted single error function and two Gaussians to mimic profiles of x-ray absorption spectra (figure 2(B)). Then, we obtained a 160 000 by 100 synthetic spectrum matrix X ₀ by multiplying W _g and ${{\boldsymbol{H}}}_{g}^{T}$ and adding Gaussian noise with SN ratio of 5.62dB. We also created a spectrum mimicking the x-ray absorption spectrum of the Li_0.5CoO₂ reference material in the same manner.

Zoom In Zoom Out Reset image size

To confirm the consistency of the results obtained with the one-sided ONMF and other methods, we used the k-means method to analyze the 2D-XAS data and compared its results with those of our method. The k-means method is a standard clustering method that classifies data into k clusters. To extract individual spatial domains with different spectral profiles, it was formulated as an optimization problem [32]:

$$\begin{eqnarray}&&\mathop{\min }\limits_{\{{{\boldsymbol{h}}}_{j}\},\{{w}_{{ij}}\}}\sum _{i=1}^{m}\sum _{j=1}^{r}{w}_{{ij}}\parallel {{\boldsymbol{x}}}_{i}-{{\boldsymbol{h}}}_{j}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 12 }$$

where x _i is the x-ray absorption spectrum at the i-th pixel, and h _j denotes the mean absorption spectrum of the j-th domain. w_ij is a one-hot representation that takes either 1 or 0 depending on whether or not x _i belongs to the j-th domain. The vector denoted by ${{\boldsymbol{w}}}_{i}^{k-{means}}$ represents the spatial structure of the i-th domain .

To quantitatively compare results obtained with the one-sided ONMF and those with the k-means method, we calculated the cross correlation of all combinations of ${{\boldsymbol{w}}}_{i^{\prime} }$ of the one-sided ONMF and ${{\boldsymbol{w}}}_{i}^{k-{means}}$,

$$\begin{eqnarray}&&{\mathrm{Corr}}_{i^{\prime} ,i}=\displaystyle \frac{{{\boldsymbol{w}}}_{i^{\prime} }^{T}{{\boldsymbol{w}}}_{i}^{k-{means}}}{\sqrt{\parallel {{\boldsymbol{w}}}_{i^{\prime} }{\parallel }_{2}\parallel {{\boldsymbol{w}}}_{i}^{k-{means}}{\parallel }_{2}}}.\end{eqnarray} \tag{ 13 }$$

To visualize the reaction inhomogeneity, Nakamura et al obtained two-dimensional maps of LC of Li_x CoO₂ in the model composite electrode from the PTE of the x-ray absorption spectrum at each pixel [9]. Here, we confirmed the consistency of the LC map and the domain patterns obtained by the one-sided ONMF. Because the LC has a one-to-one correlation with the PTE between 7728.5[eV] < PTE < 7730.2[eV], instead of comparing the LC map and the domain patterns, we compared the PTE map with the domain patterns.

First, we determined the PTE by fitting a Gaussian distribution to the x-ray absorption spectrum at each pixel and obtained two-dimensional maps of the PTE. Next, we classified pixels on the obtained map into three sets, denoted by D1, D2 and D3, which have PTEs in the ranges of (E₂, ∞), (E₂, E₁) and (0, E₁), respectively. E₁ and E₂ are thresholds for classifying all of the pixels into D1, D2, or D3.

To quantitatively compare results obtained with the one-sided ONMF and those with the PTE, we calculated the cross correlation of w _i of the one-sided ONMF and pixels of Dk using equation (13). In this comparison, the threshold E₁ was determined by maximizing the cross correlation between the pixels of D3 and the most similar w _i to the pixels of D3, while the threshold E₂ was determined by maximizing the cross correlation between the pixels of D1 and the most similar w _i to the pixels of D1.

Furthermore, to confirm the consistency between results obtained with the one-sided ONMF and with respect to the PTE, we compared the PTEs of the three domains of the one-sided ONMF that were most highly correlated with D1, D2 and D3 with the values of the thresholds E₁ and E₂.

In order to reveal the physical property changes in the respective domains, we performed Bayesian spectroscopy [26–28] on the extracted x-ray absorption spectra (XAS). Bayesian spectroscopy can decompose XAS into an absorption edge step structure, a white line (WL) characterized with the PTE, pre-edge structures with weak absorption intensities that appear on the low-energy side of the step structure, as well as and other peak structures [28]. Furthermore, by using the Bayes free energy as an information criterion [26], we can estimate the number of spectral components for the pre-edge and other peak structures through this data-driven science approach.

To decompose the XAS, equation (14) was used as a model that is a function of photon energy E.

$$\begin{eqnarray}&&{{ \mathcal F }}_{K}(E;{\boldsymbol{\theta }}):= {y}_{0}+{ \mathcal S }(E;{{\boldsymbol{\theta }}}_{{ \mathcal S }})+\sum _{k=1}^{K}{{ \mathcal G }}_{k}(E;{{\boldsymbol{\theta }}}_{k}),\end{eqnarray} \tag{ 14 }$$

where bold symbols denote the spectral parameters for the respective spectral components and y₀ is an energy-independent baseline to correct for zero. The spectral function ${ \mathcal S }(E;{{\boldsymbol{\theta }}}_{{ \mathcal S }})$ consists of a gentle step and a WL peak, as defined in equation (15).

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal S }(E;{{\boldsymbol{\theta }}}_{{ \mathcal S }}) & := & H\left[\displaystyle \frac{1}{2}+\displaystyle \frac{1}{\pi }\arctan \left(\displaystyle \frac{E-{E}_{0}}{{\rm{\Gamma }}/2}\right)\right]\\ & & +{ \mathcal V }(E;A,{E}_{\mathrm{WL}},\omega ,\eta ),\end{array}\end{eqnarray} \tag{ 15 }$$

where the first term is a step structure (an arctangent function) [13] characterized by an intensity H, the transition energy of the x-ray absorption edge E₀, and a broadening factor Γ. The WL giving the PTE is represented by the second term in equation (15), the pseudo-Voigt function [33] parameterized by an integrated intensity A, the transition energy E_WL, the spectral width ω, and the mixing ratio η between Lorentz and Gaussian shapes. Since the WL appears near the absorption edge energy E₀, a normal distribution with a standard deviation of 4.5 eV is introduced for a prior probability of the energy difference between E₀ and E_WL. The pre-edge and other peak structures are simplified and represented as a sum of Gaussian shapes ${{ \mathcal G }}_{k}(E;{{\boldsymbol{\theta }}}_{k})$ (k = 1, ⋯ ,K), as in the second term of equation (14).

In Bayesian spectroscopy [26, 27], we apply the Bayes' theorem [34] to the spectral decomposition analysis of the spectral data D := {(E₁, y₁), ⋯ ,(E_i , y_i ), ⋯ ,(E_N , y_N )}. The model ${{ \mathcal F }}_{K}({E}_{i};{\boldsymbol{\theta }})$ in equation (14) is specified by a parameter set θ , and according to the causality, the cause θ characterizes the resultant D . The joint probability of θ and D , P( θ , D ), can be expanded as P( θ , D ) = P( D ∣ θ )P( θ ) = P( θ ∣ D )P( D ), while Bayesian inference evaluates the posterior (conditional) probability distribution of the cause θ in equation (16) under the condition D given.

$$\begin{eqnarray}&&P({\boldsymbol{\theta }}| {\boldsymbol{D}})=\displaystyle \frac{P({\boldsymbol{D}}| {\boldsymbol{\theta }})P({\boldsymbol{\theta }})}{P({\boldsymbol{D}})},\end{eqnarray} \tag{ 16 }$$

When the noise superimposed on the data ({ ⋯ , y_i , ⋯ }) is normally distributed, the likelihood term, P( D ∣ θ ), is written in the following equation:

$$\begin{eqnarray*}&&P({\boldsymbol{D}}| {\boldsymbol{\theta }})={\left(\displaystyle \frac{b}{2\pi }\right)}^{N/2}\exp \left[-{bN}{{ \mathcal E }}_{K}({\boldsymbol{\theta }})\right],\end{eqnarray*}$$

where b is a quasi-inverse temperature defined as $b:= {\sigma }_{\mathrm{noise}}^{-2}$ with the standard deviation of the imposed noise on the data, and ${ \mathcal E }({\boldsymbol{\theta }})$ is defined as follows:

$$\begin{eqnarray*}&&{{ \mathcal E }}_{K}({\boldsymbol{\theta }}):= \displaystyle \frac{1}{2N}\sum _{i=1}^{N}{\left[{y}_{i}-{{ \mathcal F }}_{K}({E}_{i};{\boldsymbol{\theta }})\right]}^{2}\end{eqnarray*}$$

and is an error function of the model ${{ \mathcal F }}_{K}({E}_{i};{\boldsymbol{\theta }})$ specified by the parameter set θ . The denominator in equation (16) is a Bayesian partition function [26], and is obtained by marginalization of the numerator terms of equation (16) in parameter space θ .

Model selection in Bayesian spectroscopy is performed using the Bayesian free energy as the information criterion [26]. In the model of equation (14), the number of components K of the Gaussian line-shape is estimated by model selection, and the Bayesian partition function becomes Z(K, b), a function of the number of components K and the quasi-inverse temperature b. The Bayesian free energy is defined as $F(K,b):= -\mathrm{ln}Z(K,b)$ [26], and the model selection [26] and the noise inference [35] can be performed as follows:

$$\begin{eqnarray*}&&(\hat{K},\hat{b})=\mathop{\arg \,\min }\limits_{K,b}F(K,b).\end{eqnarray*}$$

the maximum a posteriori (MAP) estimation of the parameter set θ is performed by posterior probability maximization as following:

$$\begin{eqnarray*}&&\hat{{\boldsymbol{\theta }}}=\mathop{\arg \,\max }\limits_{{\boldsymbol{\theta }}}\left\{\exp \left[-\hat{b}N{{ \mathcal E }}_{\hat{K}}({\boldsymbol{\theta }})\right]P({\boldsymbol{\theta }})\right\}.\end{eqnarray*}$$

We assessed the performance of our method on synthetic data with eight separate domains.

First, the rank of the factor matrices r was determined by bi-cross-validation. Figure 3(A) shows the bi-cross-validation error as a function of rank r for 4 ≤ r ≤ 9. The bi-cross-validation error was minimized at r = 8; thus, the estimated rank was confirmed to be identical to the number of ground-truth domains of synthetic data.

Zoom In Zoom Out Reset image size

When the rank of the factor matrices was set to eight (r = 8), the one-sided ONMF successfully extracted eight spatial domains with different spectral profiles from the synthetic data (figures 3(B) and (C)). A1–A8 in figure 3(B) show spatial domain patterns obtained by reshaping the column vectors of the factorized matrix W into 2D arrays, and A1–A8 in figure 3(C) indicate the absorption spectra corresponding to domains A1–A8. A comparison of A1–A8 and g1–g8 in figure 2 confirmed that the spatial domains obtained by our method were identical with the ground-truth domains g1–g8 in the synthetic data.

By using the one-sided ONMF in combination with removal of the reference signal, we sought to extract spatial domains with different spectral profiles from the real 2D-XAS data of the model composite electrode. Here, we used two data sets for the model electrode charged to 4.2 V at current rates of 9 mA cm⁻² (figure 1 (B)) and 12 mA cm⁻² (figure 1 (A)) measured by Nakamura et al [9]. To reduce the calculation time, we clipped data from the whole images, as indicated by the squares in figures 1 (A) and (B).

First, the rank of the factor matrices r was determined by bi-cross-validation. Figures 4(A) and 5(A) show the bi-cross-validation error as a function of rank r for 7 ≤ r ≤ 11. The bi-cross-validation errors were minimized at r = 10 on both data.

Zoom In Zoom Out Reset image size

Next, we examined how our method decomposed the 2D-XAS data into spatial domains and their corresponding absorption spectra when the rank of the factor matrices was set to ten (r = 10) (figures 4 and 5). B1–B10 in figures 4(B) and 5(B) show spatial domain patterns obtained by reshaping the column vectors of the factorized matrix W , and B1–B10 in figures 4(C) and 5(C) indicate absorption spectra corresponding to the ten domains B1–B10 in figures 4(B) and 5(B).

Zoom In Zoom Out Reset image size

To confirm the consistency of the results obtained with the one-sided ONMF and the other methods, the k-means method was used to classify real 2D-XAS data into several clusters. Here, we used the data from the model electrode charged to 4.2 V at a current rate of 12 mAcm⁻² (figure 1 (A)). In accordance with the results of the model selection for the one-sided ONMF, we separated the 2D-XAS data into ten clusters by using the k-means method (see figure 6).

Zoom In Zoom Out Reset image size

To quantitatively compare results obtained with the one-sided ONMF and with the k-means method, we calculated the cross correlation of the domain patterns extracted by both methods. Table 1 shows the cross correlation matrix between B1–B10 in figure 4(B) and C1–C10 in figure 6. Domain patterns B1–B8 obtained with the one-sided ONMF were highly correlated with C1–C8 obtained with the k-means method.

Nakamura et al obtained two-dimensional maps of LC of Li_x CoO₂ in the model composite electrode from the PTE of the x-ray absorption spectrum at each pixel [9]. Here, we confirm the consistency of the LC map and the domain patterns obtained by the one-sided ONMF. Because the LC can be related to the PTE by a one-to-one function, instead of comparing the LC map and the domain patterns obtained by the one-sided ONMF, we decided to compare the PTE map with the domain patterns.

Figure 7(A) shows a two-dimensional map of the PTE obtained from the data on the electrode model charged to 4.2 V at a current rate of 12 mAcm⁻² in figure 1(A). As described in section 2.6, three sets of pixels, D1, D2 and D3, into which the pixels were classified with respect to the PTE by thresholds E₁ and E₂, are plotted in figure 7(B). B8, B5, and B3 of figure 7(B) are the domains of the one-sided ONMF that are most highly correlated with D1, D2 and D3. Here, the thresholds maximizing the correlations between the domain patterns of the two methods are E₁ = 7730.18[eV] and E₂ = 7729.79[eV]. Table 2 is the cross correlation matrix of B1–B10 in figure 4(B) and D1–D3 in figure 7(B). The correlation coefficients between B8 and D1, between B5 and D2, and between B3 and D3 were 0.61, 0.34 and 0.59, respectively.

Zoom In Zoom Out Reset image size

Furthermore, to confirm the consistency between domains obtained by the one-sided ONMF and by thresholding the PTE map, we compared the PTEs of the three domains B8, B5 and B3 with the values of thresholds E₁ and E₂. The PTEs of the spectra of B1–B10 in figure 4 are listed in table 3. The PTE of domain B8, 7730.38[eV], was higher than E₁ = 7730.18[eV], while the PTE of domain B3, 7729.65[eV], was lower than E₂ = 7729.79[eV]. Moreover, the PTE of domain B5, 7729.79[eV], was equal to E₂. Thus, the PTEs of the domains extracted by the one-sided ONMF were correlated with the thresholds E₁ and E₂ for separating the PTE map into the lower, middle and higher sides. Note that the thresholds E₁ and E₂ exist in the range where the LC of Li_x CoO₂ and the PTE satisfy the one-to-one relation.

Figure 7 presents three spatial domain patterns, B3, B5, and B8, which correspond to the spatial patterns D1-3 classified by PTE. We used Bayesian spectroscopy to decompose the XAS extracted from B3, B5, and B8; the results are presented in figure 8. Colored-area spectra denote the extracted XAS by the one-sided ONMF. Bayesian spectroscopy uses the MAP values, which maximize the total posterior probability. The spectral component of ${ \mathcal S }(E;{{\boldsymbol{\theta }}}_{{ \mathcal S }})$ and the peak structures of ${{ \mathcal G }}_{k}(E;{{\boldsymbol{\theta }}}_{k})$ in equation (14) are depicted as dashed and thin solid curves, respectively. By using Bayes free energy minimization, the peaks (K in equation (14)) are estimated to number 14, 12, and 14 for the XAS extracted from domains B3, B5, and B8, respectively.

Zoom In Zoom Out Reset image size

The regressed spectra ${{ \mathcal F }}_{\hat{K}}(E;{\boldsymbol{\theta }})$, which are the sums of these spectral components, are presented by thick solid curves, in which the baseline y₀ was estimated to be small enough. The root-mean-square deviations (RMSDs) of the ${{ \mathcal F }}_{\hat{K}}(E;{\boldsymbol{\theta }})$ are in the range of 1.1 − 2.0 × 10⁻³ on the normalized absorption intensity scale, and the regressive spectra explained the target XAS well because of the sufficiently small RMSD.

Such fine spectral decomposition provides important information for understanding the charging process and its spatial inhomogeneity in lithium-ion batteries, because it is considered that the transition energy, intensity, spectral width, and their changes of each decomposed spectral component reflect the physical properties of the electronic states involved in ionic conduction. Indeed, as seen in the insets of figure 8, one can find spectral changes in the decomposed components for pre-edge region on the lower energy side of the WL, in which the ordinate is displayed on a five-times enlarged scale. Here, there are differences in the decomposed spectral components ${{ \mathcal G }}_{k}(E;{{\boldsymbol{\theta }}}_{k})$ (k = 1 ∼ 3), which will be discussed in section 4.3.

The automatic domain extraction method that we developed can identify and separate very small differences in x-ray absorption spectra. We removed the reference signal from 2D-XAS data of the model composite electrode by subtracting the x-ray absorption spectrum of Li_0.5CoO₂ reference material from the 2D-XAS data; then we applied the one-sided ONMF to the difference spectra without the reference signal.

We assessed the performance of our proposed method on synthetic data that mimic profiles of x-ray absorption spectra. Bi-cross validation was used to show that number of the ground-truth domains was correctly identified and the spatial domains obtained by our method were identical with the ground-truth domains in the synthetic data (see figures 2 and 3). Next, we assessed the performance of our method on two different real 2D-XAS datasets (see figures 1(A) and (B)). We confirmed that our method could determine the rank of the factor matrices r and extract spatial domains with different spectral profiles from these data (see figures 4 and 5). Then, we quantitatively compared the results obtained with our method and those with the k-means method on the real 2D-XAS data shown in figures 1(A). We confirmed that almost all of the domain patterns obtained with our method were highly correlated with the results obtained with the k-means method (see figures 4(B) and 6 and table 1); thus, the results of the two methods were consistent. Moreover, we quantitatively compared the results obtained by our method and those by thresholding of the PTE map shown in figures 1(A). We found that the spatial patterns and the PTEs of the domains extracted by the one-sided ONMF were correlated with those of the domains obtained by thresholding the PTE map (see figure 7 and tables 2 and 3). Note that the consistency with the results of the k-means method and the PTE map had also been confirmed using other data clipped from a non-active region in figure 1(A) [36]. Finally, we applied Bayesian spectroscopy to spectra of different domains obtained by our method and decomposed the spectra of the individual domains into spectral components. We found that the spectral features of some of the decomposed spectral components differ in the individual domains (see figure 8).

These results lead us to conclude that the one-sided ONMF in combination with removal of the reference signal can extract spatial domains with different spectral profiles from 2D-XAS data showing very little variation from one position to another. Furthermore, we conclude that the individual domains and their spectra obtained by our method represent differences in physically interpretable features, including LC. The physical interpretation of the decomposed spectral components will be discussed below.

As shown in figure 6 and table 1, we demonstrated that the k-means method can also automatically extract spatial domains with different profiles of absorption spectra from the 2D-XAS data of the model composite electrode. Our method has the following advantages over the k-means method:

• 1.  
  The classification achieved by the k-means method is sensitive to the magnitude of each sample, because samples are partitioned into several clusters in which each sample belongs to the cluster with the nearest cluster center. On the other hand, the classification achieved by the NMF is invariant to the magnitude of each sample, because through the matrix factorization the magnitude of each sample is incorporated into the factor matrix W and samples are separated into several clusters in which each sample belongs to the cluster with the most directionally similar column vector in the factor matrix H. Thus, our method can achieve magnitude-invariant classification for the spectra.
• 2.  
  The domain patterns obtained with the k-means method with one-hot labels taking either 1 or 0 are not allowed to overlap others, whereas the domain patterns obtained with our method are gently allowed to overlap others because of the column orthogonality imposed on W by using the Lagrange multiplier method (equation (8)). Thus, our method can extract spatially overlapped domains.

As described in section 3.5, the XAS extracted from the spatial domains of B3, B5, and B8 were decomposed by Bayesian spectroscopy, and differences were found in some of the spectral components. In order to understand the physical meaning of these differences, the decomposed spectral components are summarized again in figures 9(a) and (b).

Zoom In Zoom Out Reset image size

In figure 7 of section 3.4, we have shown the domain patterns B8, B5, and B3 corresponding to the two-dimensional maps D1−3 with different PTEs. This change in PTE can be confirmed as a change in the transition energy of WL, which gives the maximum absorption intensity in XAS through the Bayesian spectroscopic analysis described in section 3.5.

The dashed curves in figure 9(a) show the spectral components ${ \mathcal S }(E;{{\boldsymbol{\theta }}}_{{ \mathcal S }})$ including a gentle step and an associated WL peak, as defined by equation (15): one can see that their PTEs are different. The vertical dashed lines in figure 9(a) are the MAP values of the WL-transition energy E_WL. The posterior probability distributions of E_WL can be evaluated by using Bayesian spectroscopy on each of the B3, B5, and B8 spatial domains. Their distributions are sharp, and the overlaps among them are not so obvious. Consequently, the change in E_WL between the different spatial domains can be confirmed statistically, and this result supports the result in section 3.4.

The solid curves in figure 9(a) are the sum spectra of the high-energy-side components defined by equation (17), which include ${ \mathcal S }(E;{{\boldsymbol{\theta }}}_{{ \mathcal S }})$ and peak structures ($4\leqslant k\leqslant \hat{K}$) other than pre-edge components (k = 1 ∼ 3).

$$\begin{eqnarray}&&{{ \mathcal F }}_{{\rm{a}}}(E;{{\boldsymbol{\theta }}}_{{\rm{a}}}):= { \mathcal S }(E;{{\boldsymbol{\theta }}}_{{ \mathcal S }})+\sum _{k=4}^{K}{{ \mathcal G }}_{k}(E;{{\boldsymbol{\theta }}}_{k}),\end{eqnarray} \tag{ 17 }$$

where the indexes of ${{ \mathcal G }}_{k}(E;{{\boldsymbol{\theta }}}_{k})$ are counted from the low-energy side. Although a noticeable difference appears in the WL, these curves, ${{ \mathcal F }}_{{\rm{a}}}(E;{{\boldsymbol{\theta }}}_{{\rm{a}}})$, appear to have no noticeable differences in the low energy region of E_WL.

However, we can find spatial inhomogeneities in the electronic states considered to be relating to the ionic conduction from careful examination of the spectral components of the pre-edge structure decomposed by Bayesian spectroscopy. Figure 9(b) shows the decomposed spectral components for the pre-edge structures, which appear on the low energy side of ${{ \mathcal F }}_{{\rm{a}}}(E;{{\boldsymbol{\theta }}}_{{\rm{a}}})$.

In Li_x CoO₂, Co atoms are located in the crystal field of the octahedral coordination of oxygen atoms [37], and the peaks depicted by the solid curves, which appear at the lowest energy in figure 9(b), are considered to be a pre-edge structure due to the transition into the Co-3d t_2g state [37–39]. However, this transition does not show a significant change in transition energy with respect to the spatial domains B3, B5 and B8.

On the other hand, one can find pronounced changes in the two components decomposed at the high energy side of Co-3d t_2g as shown by the dashed and dotted curves in figure 9(b).

It has been reported that the e_g state changes depending on the state of charge [37]. However, it has also been reported that there is a hybridization between the Co-3d e_g state and the 2p orbital of oxygen [39], and such a difference in hybridization is thought to appear as a change in these two spectral components. The above implies that the electronic-state changes were detected in different spatial domains through the one-sided ONMF of the 2D-XAS data measured on the model composite electrode of Li ion batteries and spectrum decomposed by Bayesian spectroscopy. The method presented here will help us to gain a better understanding of the physical factors governing electrochemical reaction inhomogeneity.