Anal. Chem. 2003, 75, 4236-4243
Advantages of Soft versus Hard Constraints in Self-Modeling Curve Resolution Problems. Alternating Least Squares with Penalty Functions Paul J. Gemperline* and Eric Cash
Department of Chemistry, East Carolina University, Greenville, North Carolina 27858
A new algorithm for self-modeling curve resolution (SMCR) that yields improved results by incorporating soft constraints is described. The method uses least squares penalty functions to implement constraints in an alternating least squares algorithm, including nonnegativity, unimodality, equality, and closure constraints. By using least squares penalty functions, soft constraints are formulated rather than hard constraints. Significant benefits are obtained using soft constraints, especially in the form of fewer distortions due to noise in resolved profiles. Soft equality constraints can also be used to introduce incomplete or partial reference information into SMCR solutions. Four different examples demonstrating application of the new method are presented, including resolution of overlapped HPLC-DAD peaks, flow injection analysis data, and batch reaction data measured by UV/visible and near-infrared spectroscopy (NIR). Each example was selected to show one aspect of the significant advantages of soft constraints over traditionally used hard constraints. Incomplete or partial reference information into self-modeling curve resolution models is described. The method offers a substantial improvement in the ability to resolve time-dependent concentration profiles from mixture spectra recorded as a function of time. Self-modeling curve resolution (SMCR) describes a set of mathematical tools for estimating pure component spectra and composition profiles from data matrixes of mixture spectra recorded from evolving systems. Evolving systems are complex chemical systems that change in systematic, nonrandom ways as a function of time, pH, temperature, etc.1 Examples of such systems include partially resolved chromatographic peaks eluting as a function of time, flow injection analysis mixtures, industrial reaction mixtures, and protein melting curves. In each example cited, changes in one measurement order (time, pH, temperature, etc.) modulate the response of a detector, for example, a spectrometer, in a second measurement order. The result is an I × J data matrix, A, with I mixture spectra in rows measured at J wavelengths that follow the so-called “bilinear model”
A ) CPT + E
(1)
where C is an I × K matrix of K pure-component concentration * Corresponding author. (1) Diewok, J.; de Juan, A.; Maeder, M.; Tauler, R.; Lendl, B. Anal. Chem. 2003, 75, 641-647.
4236 Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
profiles in columns, and PT is a K × J matrix of pure component spectra. Measurement errors, E, are assumed to be independent and have constant variance. In these kinds of evolving systems, one or more components are typically unknown and cannot be separated or isolated from the mixture in a pure state. Examples include short-lived reactive intermediates in chemical reactions, intermediate species in equilibrium mixtures, etc. In such cases, mathematical resolution of the mixture spectra is the only way pure component profiles and spectra of these species can be obtained. The SMCR Research Hypothesis. In typical SMCR applications, there is an implied research hypothesis that is not always recognized or clearly stated by authors and practitioners: “There exists an unconstrained bilinear model with unimodal, nonnegative pure component concentration profiles and pure component nonnegative spectra that fits the data matrix of measurements obtained from the evolving system.” This hypothesis is tested by iteratively fitting a constrained model to the data set until convergence is obtained. If the hypothesis is true, the SMCR calculation will yield a final model with no active constraints. This important consequence is not always recognized by practitioners of SMCR. In practice, deviations from ideal behavior due to noise, nonideal chemical response and nonideal spectroscopic response may frequently yield final models that require active constraints. Active constraints may improve the model interpretability at the expense of an increase in the model lack-of-fit. In some cases, the lack-of-fit can be so severe that it calls into question the validity of the original hypothesis. Fitting SMCR models usually begins with selection of some initial starting point followed by iterative refinement with constrained alternating least squares steps until convergence to a stationary solution. The following basic steps are required in order to fit a constrained bilinear SMCR model to a set of measurements. 1. An initial solution or “guess” is selected by any number of popular methods. If good information is available about the concentrations profiles, an initial estimate of concentration profiles, C0, may be used. Some examples include the “needle search” method2,3 or evolving factor analysis (EFA),4-6 or alternatively, spectra,7 P0, or variables8,9 from the measurements themselves, (2) Gemperline, P. J. J. Chem. Inf. Comput. Sci. 1984, 24, 206-212. (3) Gemperline, P. J. Anal. Chem. 1986, 58, 2656-2663. (4) Maeder, M.; Zuberbuehler, A. D. Anal. Chim. Acta. 1986, 181, 287-291. (5) Gampp, H.; Maeder, M.; Meyer, C. J.; Zuberbuehler, A. D. Comments Inorg. Chem. 1987, 6, 41-60. (6) Maeder, M. Chemom. Intell. Lab. Syst. 1987, 12, 527-530. 10.1021/ac034301d CCC: $25.00
© 2003 American Chemical Society Published on Web 07/22/2003
which are hypothesized to be good approximations of pure component spectra or pure variables. 2. The initial starting point, either C0 or P0, seldom obeys the constraints imposed by the research hypothesis; thus, constrained alternating least squares steps are used to fit the initial unsconstrained solution producing better “constrained” estimates. a. Given some initial or intermediate estimate of C, find P such that P minimizes ||A - CPT|| subject to constraints on P such as P > 0, etc. b. Given some initial or intermediate estimate of P, find C such that C minimizes ||A - CPT|| subject to constraints on C such as C > 0, etc. 3. The purpose of the constraints during fitting is to ensure that the original starting solution converges smoothly and monotonically to the desired result. Several recently published papers describe in detail the mechanism for solving these types of constrained problems in a least squares manner then ensures monotonic convergence, the most widely used being the nonnegative least squares (NNLS) method of Lawson and Hansen10 and adapted by Bro and de Jong11,12 for solving SMCR problems. 4. Each application of the two alternating least squares steps (2a and 2b) produces a better estimate of the constrained concentration profiles and pure component spectra; thus, a simple iterative refinement process is used whereby these two alternating steps are repeated until no further improvement in the estimates Cn and Pn is observed or the maximum number of iterations, n, is reached. Under most conditions, nonunique SMCR solutions are obtained.13 Typically, the spectroscopic techniques employed do not give unique response for each species present in the mixture, there are insufficient regions in time where components are resolved, and the constraints do not add sufficient information to completely resolve the problem. A thorough discussion of conditions that are required for unique solutions can be found in a paper by Smilde et al.14 In the absence of conditions that may lead to unique solutions, a range of feasible solutions that obey the constraints is obtained.15 Several papers have been published demonstrating methods for computation of these ranges.16,17 Hard versus Soft Constraints. In this paper, we present a new SMCR algorithm that employs “soft” constraints and compare the results to those obtained using “hard” constraints. “Hard” constraints are constraints that are strictly enforced, whereas “soft” constraints allow small deviations from constrained values. Most, if not all, existing SMCR algorithms in use to date employ hard constraints. For many data sets, we have observed that the use (7) Sanchez, F. C.; Toft, J.; van den Bogaert, B.; Massart, D. L. Anal. Chem. 1996, 68, 79-85. (8) Schostack, K. J.; Malinowski, E. R. Chemom. Intell. Lab. Syst. 1989, 6, 2129. (9) Windig, W.; Heckler, C. E.; Agblevor, F. A.; Evans, R. J. Chemom. Intell. Lab. Syst. 1992, 14, 195-207. (10) Lawson, C. L.; Hanson, R. J. Solving Least Squares Problems; Prentice Hall: Englewood Cliffs, NJ, 1974. (11) Bro, R.; De Jong, S. J. Chemom. 1997, 11, 393-401. (12) Bro, R.; Sidiropoulos, N. D. J. Chemom. 1998, 12, 223-247. (13) Tauler, R.; Kowalski, B. J. Chemom. 1995, 9, 31-58. (14) Smilde, A. K.; Hoefsloot, H. C. J.; Kiers, H. A. L.; Bijlsma, S.; Boelens, H. F. M. J. Chemom. 2001, 15, 405-411. (15) Gemperline, P. J. Anal. Chem. 1999, 71, 5398-5404. (16) Borgen, O. S.; Kowalski, B. R. Anal. Chim. Acta. 1985, 174, 1-26. (17) Borgen, O. S.; Davidsen, N.; Zhu, M.; Oeyen, O. Mikrochim. Acta 1986, 2, 63-73.
Figure 1. Schematic illustration of the row-wise fitting algorithm for finding rows of P.
Figure 2. Schematic illustration of the row-wise fitting algorithm for finding rows of C.
of hard constraints leads to final solutions with many active constraints. These types of solutions may sometimes exhibit distorted composition profiles and spectra. For example, with soft constraints, small negative deviations are allowed when nonnegative constraints are employed, thereby minimizing the number of active constraints and the impact of noise and nonideal response. Compared to hard constraints, SMCR models with soft constraints often have fewer active constraints, hence, lower model lack-of-fit, and are more likely to fit the original research hypothesis. There exists an unconstrained bilinear model with unimodal, nonnegative pure component concentration profiles and pure component nonnegative spectra that fits the data matrix of measurements. The new algorithm uses least squares penalty functions to implement constraints during alternating least squares steps; hence, it is called P-ALS. In the remaining sections of this paper, we describe the P-ALS algorithm in detail, including how penalty functions can be used to construct equality constraints,18 nonnegativity constraints, unimodality constraints, and closure constraints18 using the “row-wise” ALS method of Bro.19 After presenting the details of the new algorithm, we present four compelling examples using simulated and real data to illustrate the significant advantages of soft constraints over hard constraints in SMCR algorithms. The P-ALS Algorithm. The P-ALS algorithm uses the same four steps outlined above to alternatively estimate pure component spectra, P, and concentration profiles, C; however, the two alternating problems are solved in a row-wise fashion with least squares penalty functions to implement constraints. Bro showed that finding the optimum least squares estimate of individual rows in P (row-wise estimation of P or column-wise estimation of PT) can be performed independently for rows in the same mode.19 Thus, steps 2a and 2b outlined above become: a. Given some initial or intermediate estimate of C, for every j, find p:,j such that p:,j minimizes ||a:,j - CpT:,j|| subject to constraints on p:,j such as p:,j ) g:,j. The global solution to P is obtained by solving the row-wise subproblem j times. The procedure is illustrated schematically in Figure 1. By transposing the problem, the same algorithm can be used to solve the row-wise estimation of CT, as shown in Figure 2. (18) Van Benthem, M. H.; Keenan, M. R.; Haaland, D. M. J. Chemom. 2002, 16, 613-622. (19) Bro, R. University of Amsterdam, 1998.
Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
4237
b. Given some initial or intermediate estimate of PT, for every i, find ci,: such that ci,: minimizes ||aTi,: - Pci,:T|| subject to constraints on ci,:, such as ci,: ) gi,:. Equality Constraints. Approximate equality constraints can be implemented for least squares problems y ) Xb via penalty functions.18,19 Notice that both row-wise P-ALS subproblems (a) and (b) described can be generalized to this form. Suppose the equality constraints bi ) gi are desired, where the elements gi represent the desired goals for selected coefficients bi of b. By augmenting y with g and X with a matrix, H, of appropriately positioned ones and zeros, one or more coefficients of the normal least squares solution vector, b, can be forced to conform to the desired values, gi (see numerical examples under the section Nonnegativity Constraints).
[ ] [ ]
y X = b λg λH
(3)
The penalty function weighting factor, λ, can be adjusted to small values when soft constraints are desired, or it can be set very large to give hard constraints. To properly weight problems of different sizes and different measurement scales, the penalty function weighting factor, λ, can be adjusted relative to the norm of X; for example, 0.01 × norm (X) for soft constraints or 10 × norm (X) for hard constraints. If one or more pure component spectra are known a priori, this information can be included as equality constraints on PT in step (a) of the P-ALS method.18 If reference data is available for concentrations, say for example, aliquots of a reaction mixture are analyzed at selected time intervals by HPLC, these sparsely sampled reference data can be used as equality constraints on C in step (b) of the P-ALS method. Nonnegativity Constraints. A simple modification of P-ALS steps 2a and 2b can be used to impose approximate nonnegativity constraints using least squares penalty functions and equality constraints. First, the standard unconstrained least squares problem is solved. The resulting coefficients, bi, are inspected. For any coefficients, bi, that are negative, an equality constraint, bi ) 0, is set. With a large penalty value, the resulting solution converges to the same result that would be obtained using hard constraints, such as those obtained from algorithms NNLS10 or FNNLS.11 An example illustrating the implementation of nonnegative constraints is shown below. First, the unconstrained least squares solution to a sample problem is given.
[ ][
0.9218 0.7382 y) X) 0.1763 0.4057
0.4451 0.9318 0.4660 0.4186
0.8462 0.5252 0.2026 0.6721
][
0.8381 0.2660 0.0196 bˆ 0.8163 0.6813 -0.0595 0.3795
]
Noting that coefficient b3 is negative, the nonnegatively constrained least squares problem is solved as shown below with equality constraints and penalty weighting function λ ) 10, giving the approximate nonnegatively constrained solution, bP-ALS(λ:10) ) [0.2749, 0.7645, -0.0003]. If hard constraints are imposed by use of algorithm NNLS, the solution bNNLS ) [0.2749, 0.7643, 0] is obtained. This solution compares favorably to the solution obtained with penalty weight λ ) 100 bP-ALS(λ: 100) ) [0.2749, 0.7643, -3 × 4238
Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
10-6]. If soft constraints are desired, a smaller penalty weight, such as λ ) 1.0, can be used, producing the solution bP-ALS(λ:1) ) [0.2721, 0.7809, -0.0189].
[ ][
0.9218 0.7382 y ) 0.1763 X ) 0.4057 0.0
0.4451 0.9318 0.4660 0.4186 0.0
0.8462 0.5252 0.2026 0.6721 0.0
][
0.8381 0.0196 0.2749 0.6813 bˆ ) 0.7645 0.3795 -0.0003 10
]
Closure Constraints. In some cases, such as in batch reaction studies or chemical equilibrium studies, the principle of mass balance can be invoked such that the sum of all or selected species concentration profiles should equal a constant. Closure constraints can be implemented with equality constraints in the manner shown by Van Benthem, Keenan, and Haaland18 using H ) [1, 1, ... 1] and g ) [1]. Augmenting the previous example with closure constraints with a penalty weighting function λ ) 10 gives the following results with sum(bP-ALS) ) 1.0005. The example also illustrates that many different combinations of constraints can be solved in one step.
[ ][
0.9218 0.7382 0.1763 y ) 0.4057 X ) 0.0 0.0 10
0.4451 0.9318 0.4660 0.4186 0.0 10
0.8462 0.5252 0.2026 0.6721 0.0 10
][
0.8381 0.0196 0.245 0.6813 bˆ ) 0.7645 0.3795 -0.0004 10 10
]
For some chemical systems, a weighted sum of several species would be expected to give a constant, for example, as in the reaction 2A f B f C. Here, the appropriate closure constraint would be H ) [1/2 1 1] and g ) [1]. Unimodality Constraints. Further modification of P-ALS step 2b can be used to impose approximate unimodality constraints using least squares penalty functions and equality constraints. In a fashion similar to the nonnegativity constraint implementation described above, the standard unconstrained least squares problem is solved first. The resulting concentration profiles are inspected as a function of time to find the global peak maximums, one for each profile. Searching in both the forward and reverse directions from the global peak maximums for each constituent, unimodality constraints must be added at time i + 1, gi+1,j ) ci,j, or at time i - 1, gi-1,j ) ci,j, if secondary local maximums are encountered, for example, ci,j < ci+1,j or ci,j < ci-1,j. An example illustrating the implementation of unimodality constraints is shown below. First, the unconstrained least squares solution to a sample problem is given.
[ ][
0.5534 0.2920 y) X) 0.8580 0.3358
0.1870 0.9913 0.7120 0.8714
0.4796 0.4960 0.2875 0.0609
]
[ ]
0.2625 0.0766 0.1863 bˆ ) 0.4085 0.9171 0.7832 0.1233
Suppose it is suspected that a peak maximum occurs at b2. Noting the coefficient b3 in the unconstrained solution is too large, the unimodal constrained least squares problem is solved, as shown
below, with the equality constraint b3 ) 0.4085 and penalty weighting function λ ) 10, giving the approximate unimodalconstrained solution, bP-ALS(λ:10) ) [0.1720, 0.5908, 0.4102].
[ ][
0.5534 0.2920 y ) 0.8580 X ) 0.3358 4.0849
0.1870 0.9913 0.7120 0.8714 0.0
0.4796 0.4960 0.2875 0.0609 0.0
][
0.2625 0.1720 0.1863 0.9171 bˆ ) 0.5908 0.4102 0.1233 10
]
After several iterations of the ALS algorithm, the result converges to bP-ALS(λ:10) ) [0.1408, 0.5311, 0.5323]. Software Implementation of P-ALS. A program called GUIPRO with a graphical user interface has been written to implement the above techniques in the most general way possible. The program runs under Matlab 6.0 or higher. Copies of this program may be obtained from the author by sending an email to
[email protected]. The user interacts with the program using the mouse to click buttons to perform self-modeling curve resolution (SMCR). The program can load spectroscopic data from Matlab files or tab-separated data from EXCEL spreadsheets. GUIPRO’s menu structure is designed to guide users through the correct sequence of steps needed to perform SMCR, including importing data, selecting subranges for analysis, removing outliers by inspecting PCA scatter plots, correcting baseline shifts, normalization of spectra, selecting number of PCA factors, performing curve resolution, and saving and reloading results at a later time. Performing curve resolution with the “const. P-ALS” option selected brings up the following graphical user interface for setting and controlling the behavior of P-ALS. Using this graphical user interface, constraints can be turned on or off separately for each constituent by selecting the appropriate check boxes. For concentration profiles, nonnegativity, unimodality, and closure constraints may be selected. For spectral profiles, nonnegativity constraints may be selected. The penalty function weighting factors for concentration profiles and spectral profiles may be set separately using the sliders shown on “global options”. When equality constraints are desired, auxiliary reference data must be loaded from the Matlab workspace. Equality constraints are automatically generated for any portion of concentration profiles where auxiliary concentration data are loaded and equality constraints are automatically generated for any portion of pure component spectra where auxiliary reference spectra are loaded. The auxiliary data may be sparsely sampled or sampled at different times, as compared to the sampling times for spectral data. When the auxiliary concentration data sample times do not coincide exactly with the sample times for spectral data, the software offers two options for synchronizing the auxiliary concentration data to sample time for spectral data. When the “discrete” button is checked, nearest neighbors between auxiliary concentration sample times and spectral sample times are used. When the “continuous” button is checked, the auxiliary concentration data are interpolated to match the spectral sample times. If pure component spectra are available for one or more spectral profiles, these may be loaded as auxiliary reference data for spectral profiles. As with auxiliary concentration data, interpolation or nearest neighbors are automatically used when the sample
Figure 3. Screen shot of P-ALS options window.
intervals of auxiliary reference spectra do not exactly match the sample intervals for the data set to be analyzed. EXPERIMENTAL SECTION Four different examples are provided to show the advantages of using soft constraints instead of hard constraints in SMCR problems. (1) A simulated HPLC-diode array data set is shown in which overlapped peaks are mathematically resolved at different noise levels. This data set is used to illustrate problems caused by noise and how the use of soft constraints significantly improves the quality of the results. (2) A flow injection analysis/diode array data set is analyzed showing the resolution of Bi3+, BiCl2+, BiCl2+, BiCl3, BiCl4-, BiCl52-, and BiCl63- in a bolus of Bi(ClO4)3 injected into a flowing stream of HCl.20 This data set is used to show the advantages of using reference spectra and soft equality constraints to improve the results of SMCR solutions. (3) The characterization UV/visible ATR spectra acquired in situ during a batch reaction of salicylic acid with acetic anhydride.21 This data set is used to illustrate the use of temperature data and soft equality constraints to identify a temperature-dependent shift in a UV/visible absorption band. (4) The characterization of NIR spectra acquired in situ during the batch reaction of butanol and acetic anhydride. This data set is used to illustrate the use of sparsely sampled concentration data to improve the quantitative results from SMCR analysis of batch reaction data. Simulated HPLC/DA Data. The chromatograms and spectra shown in Figure 4 were used to generate simulated threecomponent mixtures with overlapping Gaussian chromatography peaks to illustrate the effects of measurement error on the quantitative accuracy of solutions produced by SMCR with hard and soft constraints. The data matrix was scaled to have a maximum signal intensity of 1.0. Random measurement errors (20) Gemperline, P. J.; Hamilton, J. C. J. Chemom. 1989, 3, 455-461. (21) Ma, B.; Gemperline, P. J.; Cash, E.; Bosserman, M.; Comas, E. J. Chemom., 2003, in press.
Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
4239
Figure 4. Overlapped chromatographic peaks and pure component spectra used to generate simulated HPLC-DAD data set.
were simulated by adding normally distributed random numbers with a mean of zero and a standard deviation of 0.02 to the simulated data matrix. Flow Injection Analysis Data. In acidic solutions, bismuth exists in equilibrium with chloride22 in the complexes Bi3+, BiCl2+, BiCl2+, BiCl3, BiCl4-, BiCl52-, and BiCl63-. A previously published flow injection analysis/diode array data set was obtained by injecting a bolus of Bi(ClO4)3 into a flowing stream of HCl20. Using evolving factor analysis, the composition profiles and pure component spectra of the various components were estimated. The apparatus was designed to control the amount of mixing before the sample bolus reached the detector so that at the peak apex, the concentration of the carrier stream was reduced to nearly zero, the concentration of Bi3+ was at a maximum, and minimal complexation with Cl- was observed. In the leading edge and tail of the bolus, the concentration of Bi3+ was reduced to nearly zero, and the concentration of Cl- was at its maximum, so complete complexation to BiCl63- was observed; thus, a data set spanning the full range of all seven bismuth species was observed between the leading edge of the bolus, the apex, and in the trailing edge of the bolus. Full experimental details are available in a previous publication.20 UV/Visible Batch Reaction Data. UV/visible spectra were measured during the acetylation of salicylic acid (SA) with acetic anhydride (AA) with phosphoric acid as a catalyst in an AutoMate automated reaction calorimeter (Hazard Evaluation Laboratory Limited, Hertfordshire, England) controlled by WinISO software. A 2.5-mL portion of o-phosphoric acid (Fisher Scientific, Fair lawn, NJ) was added to a 500-mL container of acetic anhydride (Fisher Scientific, Fair lawn, NJ) and mixed completely. A 30-mL portion of this mixture was used as the initial reactor charge in the 50mL reactor vessel, and the reactor system was equilibrated for 40 min at a temperature of 65 °C. A 7.065-g aliquot of salicylic acid (Fisher Scientific, Fair lawn, NJ) was dosed into the glass reactor vessel after equilibration, and the mixture was allowed to react for 1 h. The reactor jacket temperature was held constant with an oil circulator system at 55 °C. The control software was programmed to hold the reactor temperature at 65 °C; however, a large endotherm from the dissolution of SA followed by a large exotherm produced by the reaction of SA with AA to form acetylsalicylic acid (ASA) caused temperature fluctuation of 16.4 (22) Spivakov, B. Y.; Stoyanov, E. S.; Gribov, L. A.; Zolotov, Y. A. J. Inorg. Nucl. Chem. 1979, 41, 453-455.
4240 Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
Figure 5. Resolved composition profiles for three overlapped peaks using (A) soft constraints and (B) hard constraints: estimated profiles (s), actual profiles (- - -).
°C. Spectra were recorded every 30 s by a four-channel fiber optic CCD UV/visible spectrograph from Equitech International, a 1024 × 512 element thermoelectrically cooled CCD array, and a Xe flash lamp. The spectrograph was equipped with an ATR fiberoptic probe having a hemispherical sapphire three-bounce ATR crystal with an effective total path length of ∼1.5 microns. NIR Batch Reaction Data. NIR spectra were measured during the acetylation reaction of butanol with acetic anhydride in a 250-mL jacketed round-bottom flask. The temperature of the reaction was thermostated with a recirculating heater/chiller at 50 °C. The reactor was initially charged with 30.0 mL of butanol (Fisher Scientific, Fair lawn, NJ) and allowed to equilibrate for ∼30 min. A 31.0-mL portion of acetic anhydride (Fisher Scientific, Fair lawn, NJ) was added to the reactor along with 2 mL of concentrated sulfuric acid (Fisher Scientific, Fair lawn, NJ) as a catalyst. Spectra were recorded every 30 s (average of 10) from 1100 to 2500 nm using a FOSS-NIRSystems 6500 dispersive spectrometer and an Optiprobe fiber-optic probe operated in the transflectance mode. The path length was adjusted to 0.5 mm, giving an effective total path length of 1.0 mm. RESULTS Simulated HPLC-DA Data. The results obtained using hard versus soft constraints for the simulated HPLC-DA data are shown in Figure 5. Except for the magnitude of the penalty weight used, identical processing conditions, including nonnegativity and unimodality constraints, were used for both examples. For soft constraints, a penalty weight of 1.0 was used. For harder constraints, a penalty weight of 20 was used. Only one example is shown here, rather than a broad range of simulated conditions, because it has been well-established by others that SMCR can give quantitative results under a wide range of peak heights and resolution.2,3,8,9 Because soft constraints were used in Figure 5A, small negative and positive deviations were allowed in the resolved profiles. In the presence of harder constraints, these small negative excursions were not permitted, as shown in Figure 5B. Harder constraints eliminated these small negative excursions in the tails of the resolved peaks, causing the peak estimates to be too narrow and slightly distorted. The distortion observed in the resolved profiles and resolved spectra significantly impact the quantitative accuracy of the integrated signal for the three peaks. Integrated peak areas for the each of the k components were computed by
Figure 6. Estimated concentration profiles (A) and pure component spectra (B) of Bi3+ (o), BiCl2+ (0), BiCl2+ ()), BiCl3 (4), BiCl4- (3), BiCl52(f), and BiCl63- (×). Table 1. Resolved Peak Areas for Three Overlapped Peaks Using Soft and Hard Constraints peak 1
peak 2
peak 3
area % error area % error area % error av % Eerror actual soft constraints hard constraints FNNLS
150.8 153.8 182.3 171.4
1.9 20.8 13.6
184.4 79.3 -2.8 182.5 -1.0 164.7 -10.7
185.3 187.1 1.0 155.8 -15.9 187.8 1.3
1.9 12.6 8.6
summing the total signal of the resolved component formed by the outer product ckpkT (see Table 1). Errors in the estimated pure component spectra (not shown) also contribute to the peak area errors shown in Table 1. With soft constraints, the average peak signal area errors were less ∼2%, whereas the errors with hard constraints were significantly larger, >12%. When algorithm FNNLS was used to solve the same SMCR problem, (results not shown) similar distortions were observed, with average peak area error of 8.6%. Flow Injection Analysis Data. Application of SMCR to the bismuth FIA data presented several challenges. Since all seven bismuth species are encountered in the leading half of the peak and the same seven species are encountered in the trailing half of the peak, unimodal constraints cannot be applied unless the data set is split in half and the two halves areanalyzed separately. Additionally, the shape of the bolus as it elutes through the detector is not symmetrical. Instead, it exhibits a significant amount of tailing due to slower laminar flow at the walls of the tubing. As a result, significant differences were observed in SMCR results for the front half of the bolus compared to the second half of the bolus.20,23 In both halves of the bolus, principal component analysis clearly revealed the presence of seven components. In the front half of the bolus, P-ALS with soft nonnegativity and unimodality constraints clearly resolved all seven components; however, in the second half of the peak, resolution of all seven different components was incomplete as a result of poorer resolution of components BiCl52+ and BiCl63+ in the tail of the bolus. To resolve the components in the second half of the peak, the pure component spectra estimated by P-ALS analysis of the front half of the peak were used as reference spectra in soft equality constraints (λ ) 0.20) for resolution of the second half (23) Malinowski, E. R. J. Chemom. 1992, 6, 29-40.
of the peak. Using this strategy, the results for the second half of the peak converged rapidly and smoothly, allowing for complete estimation of all seven components. The results of the two P-ALS models were spliced together into one plot, shown in Figure 6A. The pure component spectra in Figure 6B compare favorably to previously published pure component spectra estimated by means of chemical equilibrium calculations.24 With soft constraints, 99.98% of the total sum of squares was modeled, and small negative deviations were permitted in the solution between 240 and 300 s. With hard constraints (results not shown, λ ) 20.0), there were more active constraints, a smaller fraction, 98.86% of the total sum of squares was modeled, and the concentration profiles for BiCl52+ and BiCl63+ in the range from 240 to 370 s were noisier. UV/Visible Batch Reaction Data. Analysis of the batch reaction spectra of the acetylation of SA is problematic because there is a large temperature fluctuation (16.4 °C) during the reaction caused by the initial endothermic dissolution of SA, followed by the exothermic reaction of SA with AA. In very dilute solutions, such temperature changes might not cause significant problems; however, this batch was performed at high concentrations near the saturation limit of the product, ASA, which is more relevant for real industrial processes. At high concentrations, small perturbations, such as a change in temperature, can cause a significant fluctuation in spectroscopic response. Principal component analysis of the data set indicates the presence of three factors, suggesting three unique spectroscopically active species are present in this reaction, whereas only two unique species are expected on the basis of known chemical and spectroscopic behavior. SMCR analysis using three factors and the usual nonnegativity constraints failed to give meaningful results. For example, the needle search shows only two local minimums and fails to show the location of three local minimums; thus, some guesswork is needed to estimate the initial location of the third peak. Regardless of the location chosen for the three local minimums, SMCR convergence was slow with a significant number of active nonnegativity constraints at the final solution, suggesting the hypothesized nonnegative model is invalid for the experimental measurements at hand; for example, strict nonnegativity constraints are not valid for this experiment. (24) Kankare, J. J. Anal. Chem., 1970, 42, 1322-1326.
Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
4241
Figure 7. Estimated pure component profiles (A) and spectra (B) for the acetylation of SA: SA (-‚-), ASA (- -), and scaled temperature profile and temperature-dependent absorption band shift (s).
Noting that there was a significant temperature fluctuation over the course of the batch, we used P-ALS to test an alternative hypothesis: there exists a temperature-dependent response in the in-situ reaction spectra. To test this hypothesis, the measured temperature profile from the batch reaction was used as auxiliary reference data to set up soft equality constraints for one of the three “components” detected by PCA. All other constraints (nonnegativity, unimodality, and closure) were turned off for this component. Using these settings, P-ALS converged smoothly, giving the results shown in Figure 7. Two of the profiles show clearly the disappearance of SA and formation of product, ASA. Independent verification of the results can be obtained by noting that the estimated pure component spectra of SA and ASA agree well with measured pure component spectra (not shown). The third profile closely tracks the change in temperature during the course of the batch. Its corresponding spectrum shown in Figure 7B has a derivative-like feature with a zero crossing at ∼320 nm, typical of the difference one would expect to see between a shifted and unshifted absorption band at 320 nm. This behavior is consistent with what one would expect to see for a temperaturedependent peak shift. A large temperature decrease produced a shift to longer wavelengths, whereas a small temperature increase produced a small shift to shorter wavelengths. Use of soft constraints allowed a loose model for the temperature profile (Figure 7A) and a reasonable model for the temperature shift (Figure 7B), whereas use of hard constraints (not shown) reproduced the temperature profile exactly but gave a poorly defined model for the temperature-dependent absorption band shift. NIR Batch Reaction Data. SMCR analysis of the batch NIR spectra of the acetylation of butanol is problematic, because the intrinsic rank of the data matrix is three, not four, whereas there are four different chemical species present during the reaction. This SMCR problem can be solved by combining data from multiple batches;25 however, the purpose of this example is to show how inclusion of additional information, in the form of soft equality constraints using sparse reference concentration data, can significantly improve the SMCR results. In Figure 8A, component profiles estimated by SMCR with nonnegativity constraints and (25) Izquierdo-Ridorsa, A.; Saurina, J.; Hernandez-Cassou, S.; Tauler, R. Chemom. Intell. Lab. Syst. 1997, 38, 183-196.
4242 Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
Figure 8. SMCR estimated component profiles for the acetylation of butanol without reference data (A) and with equality constraints and sparse concentration data (B): starting material (s), added reagent (- - -) and product (-‚-).
unimodality constraints are shown. The profiles of the starting material, AA, and added reagent, butanol, are estimated incorrectly. This occurs because in the absence of additional information, SMCR nonnegativity constraints tend to force each profile to zero at some point in time during the reaction. In Figure 8B, the estimated spectra of the three species are (1) the spectrum of pure acetic anhydride (solid line), (2) the spectrum of a mixture of butanol and acetic anhydride (dashed line), and (3) the spectrum of a mixture of butyl acetate and acetic acid. The error occurred because the second spectrum was incorrectly estimated as a mixture of butanol and acetic anhydride instead of pure butanol. In Figure 8B, reference information was added, and soft equality constraints were formulated using the reference information and P-ALS. At t ) 0, equality constraints were to give acetic anhydride mole fraction equal to 1. At t ) 9 min, the point coinciding with the introduction of butanol, equality constraints were set to give acetic anhydride and butanol and mole fractions equal to 0.5. By setting these soft equality constraints with sparse information, the concentration profiles converged to the result shown in Figure 8B. The estimated pure component spectra of the three profiles (not shown) accurately reproduce the pure component spectra of (1) acetic anhydride, (2) butanol, and (3) the sum of the pure component spectra of butyl acetate plus acetic acid. Although there are other useful methods for solving this problem, this example shows that sparse reference information
in the concentration domain can be successfully used with soft equality constraints to correctly estimate the concentration of mixtures. CONCLUSIONS Using four different examples, a new SMCR algorithm called P-ALS has been demonstrated. By using soft constraints with penalty least squares functions, small deviations from the constraint settings are permitted. This approach is generally useful and in many circumstances may yield results that are better than those obtained with hard constraints. Examples were given showing that allowing small negative deviations for noise in nonnegatively constrained problems can significantly improve the accuracy of SMCR solutions. A software implementation has been described that is very flexible and can incorporate a variety of physically useful constraints and partial or incomplete reference information. For example, when one or more pure component spectra are available, soft equality constraints allow small deviations in estimated spectra. This may be especially useful in highly concentrated mixtures measured for example by NIR spectroscopy where spectral features are significantly modified in mixtures
compared to pure components. Most significantly, incomplete or sparse reference concentration information can be used in soft equality constraints to accurately estimate the concentration profiles. Copies of the P-ALS program for Matlab may be obtained from the author by e-mail. ACKNOWLEDGMENT This material is based upon work supported by the National Science Foundation under Grant no. 0201014. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. The Measurement and Control Engineering Center (MCEC) at the University of Tennessee, Knoxville, is also recognized for partial support of this work.
Received for review March 26, 2003. Accepted June 13, 2003. AC034301D
Analytical Chemistry, Vol. 75, No. 16, August 15, 2003
4243
