Neural Parameter Regression for Explicit Representations of PDE Solution Operators

Authors: konrad  mundinger
max  zimmer
sebastian  pokutta
Venue: ICLR_cc_2024_Workshop_AI4DiffEqtnsInSci
Published: 2024-02-02 17:55:14

Abstract: We introduce Neural Parameter Regression (NPR), a novel framework specifically developed for learning solution operators in Partial Differential Equations (PDEs). Tailored for operator learning, this approach surpasses traditional DeepONets (Lu et al., 2021a) by employing Physics-Informed Neural Network (PINN, Raissi et al., 2019) techniques to regress Neural Network (NN) parameters. By parametrizing each solution based on specific initial conditions, it effectively approximates a mapping between function spaces. Our method enhances parameter efficiency by incorporating low-rank matrices, thereby boosting computational efficiency and scalability. The framework shows remarkable adaptability to new initial and boundary conditions, allowing for rapid fine-tuning and inference, even in cases of out-of-distribution examples.

Main paper: 
 introduction: 
 Partial Differential Equations (PDEs) are central to modeling complex physical phenomena across diverse application areas. Traditional approaches often rely on numerical methods due to the scarcity of closed-form solutions for most PDEs. The emergence of Physics-Informed Neural Networks (PINNs) (Raissi et al., 2019) has revolutionized this domain. PINNs integrate the governing PDEs, along with initial and boundary conditions, into the loss function of a Neural Network (NN), enabling self-supervised learning of approximate PDE solutions.While PINNs are mostly employed to solve single instances of (possibly parametric) PDEs, Physics-Informed DeepONets (Wang et al., 2021b) have extended their capabilities to the learning of solution operators, i.e., mappings from initial conditions to solutions. We introduce Neural Parameter Regression (NPR), advancing operator learning by combining Hypernetworks (Ha et al., 2017) with Physics-Informed operator learning. Our approach significantly deviates from existing methods (de Avila Belbute-Peres et al., 2021; Zanardi et al., 2023; Cho et al., 2023) by parametrizing the output network as a low-rank model which by design satisfies the initial condition and is not required to learn an identity mapping. We refer to Section A.1 for a detailed discussion of related work.Our main contribution is a novel combination of Hypernetwork approaches with PINN techniques in the setting of operator learning. We experimentally demonstrate the ability of our approach to accurately capture (non-)linear dynamics across a wide range of initial conditions. Additionally, we showcase the remarkably efficient adaptability of our model to out-of-distribution examples. Section 2 introduces related concepts and our proposed method, Section 3 describes the training procedure and presents experimental results, and Section 4 concludes the paper. Detailed discussions and extended results can be found in the appendix. 
methods: 
 PINNs have emerged as a powerful tool for approximating solutions of PDEs by employing automatic differentiation to construct a loss function which is zero if and only if the NN solves the PDE.Hypernetwork H Φ u 0 (x 1 ) . . .u 0 (x m ) u 0 ∈ X h 1 h 2 . . . h p Target Network T θ = T HΦ(u0)T θ (t, x) t ⊕ x The DeepONet framework (Lu et al., 2021a) employs NNs to approximate operators in Banach spaces. It evaluates an operator G through G(u)(y) for functions u and vectors y, using a Branch Network G branch for encodings of u and a Trunk Network G trunk for encodings of y, combined via the dot product. Initially designed for supervised operator learning, its extension to self-supervised learning of PDE solution operators is achieved through Physics-Informed DeepONets (Wang et al., 2021b) . On the other hand, Hypernetworks (Ha et al., 2017) are a meta-learning approach, where one NN (the Hypernetwork) learns the parameters of another NN (the Target Network). We refer the reader to the appendix for more in-depth descriptions of these techniques.We consider initial boundary value problems (IBVPs) of the formEQUATIONHere, Ω ⊂ R d for some d ∈ N, T > 0 is the final time, u 0 the initial condition, B describes the boundary condition, and N is an operator typically composed of differential operators and forcing terms. Precisely, we are interested in approximating the Solution Operator, i.e., the mappingEQUATIONbetween the infinite-dimensional Banach spaces X and Y, where K is compact. In the following we assume that Equation 1 is well-posed, i.e., that the solution (in a suitable sense) exists and is unique.We propose (Physics-Informed) Neural Parameter Regression (NPR), a novel architecture for learning solution operators of PDEs. The architecture is depicted in Figure 1 . It consists of a Hypernetwork H Φ with parameters Φ, which takes the role of the Branch Network in the DeepONet. By mapping a discretization of a function u 0 to a p-dimensional vector, we compute the parameters θ of another NN denoted by T θ . Further, we pass (t, x) through this learned network and obtain an approximation of G(u 0 )(t, x), such that H Φ (u 0 ) explicitly approximates the function G(u 0 ).A solution operator of an IBVP faces the challenge of implicitly learning an identity mapping, as G(u 0 )(0, •) = u 0 must hold for all u 0 . Since learning the identity is a notoriously hard task (Hardt & Ma, 2017) , an essential part of our approach is to enforce the initial conditions in the Target Networks by parametrizing the output as the deviation from the initial condition over time. This approach is known as Physics-Constrained NNs (Lu et al., 2021b) and we refer to Section A.6 for a detailed description. Precisely, we reparametrize the output of T θ asEQUATIONConceptually, the network H Φ approximates a mapping between the function spaces X and Y.Effectively however, it is a mapping from R denc to R p , where d enc is the encoding dimension (number of sensors in Lu et al. (2021a) ). Since the input and output dimensions are low, the majority of parameters is in the hidden weights, which scale in O(n hidden d 2 hidden ), where n hidden is the number of hidden layers and d hidden is the hidden dimension. By parametrizing the hidden weights as a product of two matrices A i and B i of maximal rank r (Hu et al., 2022), we can reduce the number of parameters to O(n hidden d hidden r). More details can be found in Section A.8. 
experimental setup and results: 
 Given an IBVP (Equation 1), we define a compact set K ⊂ X of initial conditions and specify a sampling procedure. We parametrize the Hypernetwork H Φ using a Multi-Layer Perceptron (MLP). The input to H Φ is a batch of discretized initial conditions u 0 ∈ K. From the output, we compute the parameters of the Target Network T θ , itself an MLP with low-rank hidden weights. The input to T θ is a batch of space-time coordinates (t, x). As per Equation 3, the output of T θ is then reparametrized to satisfy the initial condition. Using automatic differentiation, we compute the derivatives of the output of T θ w.r.t. x and t and construct the residual loss L PDE . The boundary loss L BC is computed in a supervised manner. Finally, we build the total loss L as a weighted sum of L PDE and L BC and differentiate the parameters Φ of H Φ w.r.t. L to train the Hypernetwork using Adam (Kingma & Ba, 2015) . Periodic loss weight calculations normalize the magnitudes of the updates caused by the different loss components. Section A.9 contains a full description of the algorithm.We evaluate NPR on two one-dimensional example problems, the heat equation and Burgers equation. We perform ablations regarding the hidden dimension and the rank of the weights in the Target Network and compare to Physics-Informed DeepONets. Finally, we demonstrate the adaptability of our model by evaluating and fine-tuning on out-of-distribution initial conditions. All experiments are carried out on a single NVIDIA Tesla V100 GPU. The significantly higher training time for the heat equation is due to the presence of second order derivatives and the more involved sampling scheme.Heat equation. We consider the heat equation with constant Dirichlet boundary conditions on Ω = [0, 1] and T = 1. The set K is parametrized by Fourier polynomials with a fixed number of coefficients and bounded absolute value, see Section A.2.1 for details. We compute a reference solution using a finite difference scheme for twelve different initial conditions from the proposed distribution, and compute the L 1 , L 2 , and L ∞ errors between reference solution and model output.We fine-tune a trained model to the out-of-distribution initial condition u 0 (x) = 5x + 3 sin(4πx).To achieve this, we compute the output of the Hypernetwork H Φ (u 0 ) and use it as parameters for Target Network T θ . We discard H Φ and regard T θ as a conventional PINN, training it for 200 steps. The results are shown in Figure 3 . Without fine-tuning, the model performs poorly, but after only 200 iterations (≈ 2 seconds on a AMD Ryzen 7 PRO 6850U CPU) the model adapts to the new initial condition. See Section A.10 for details on fine-tuning. . In this setting we can guarantee that no shocks occur and we can give an explicit formula for the solution, see Section A.2.2 for details. Figure 4 shows how our model captures the induced nonlinear dynamics at different output scales. Close to the formation of a shock, the performance degrades slightly, which is expected given the increasingly steep nature of the solution in this region. mensional domains. This limitation arises from the exponential increase in the number of sensors needed to meaningfully represent functions as the dimensionality of their domain grows. Further, the considered examples do not show challenging behaviour such as sharp local features. We leave the exploration of these limitations to future work.A.1 RELATED WORKOur study extends the framework of PINNs and operator learning, a domain significantly advanced by the introduction of Physics-Informed DeepONets (Wang et al., 2021b) . These models leverage automatic differentiation to learn mappings from initial conditions to solutions, demonstrating a significant extension of PINN methodologies. Subsequent enhancements to the DeepONet architecture, including a weighting scheme based on Neural Tangent Kernel (NTK) theory (Jacot et al., 2018) Zanardi et al. (2023) have also applied HyperNetworks to Physics-Informed problems and also have used LoRa weight updates to a base model. We differ from this by inherently parametrizing the output net as a low-rank model. Cho et al. (2023) have also applied Hypernetworks to PINNs and parametrized the networks via low-rank matrices, but also in the setting of parametric PDEs with fixed initial conditions. To the best of our knowledge, this work is the first to extend the concept of HyperNetworks to the setting of learning solution operators of PDEs.Equation 1 describes a wide range of initial boundary value problems (IBVPs) for PDEs. The operator B can describe many different types of boundary conditions, e.g., Dirichlet, Neumann or Robin boundary conditions. We point out that the operator N may depend on x, t and additional parameters c, i.e. N = N (u, x, t, c), where we have omitted the dependence on x, t and c in the notation for the sake of readability. We only consider solution operators on compact subsets of X , as compactness is an essential requirement for the universal approximation theorem for operators (Chen & Chen, 1995) , on which the DeepONet framework was built.A solution to Equation 1 in the classical sense is a function u : [0, T ] × Ω → R which satisfies the PDE and the initial-and boundary conditions. In this setting, the space X is usually chosen as C k (Ω, R)1 , where k is high enough to allow for evaluations of the differential operator N . The space Y is then usually chosen as C([0, T ] × Ω, R). Many other solution concepts exist, e.g. weak solutions, mild solutions or viscosity solutions. We refer to the classical literature on PDEs (Evans, 1998) for an overview of solution concepts. target network) and between 15 and 22 minutes for the Burgers equation. The difference is due to higher order derivatives in the heat equation.For the DeepONets in Table 1 we used a Branch Network with four hidden layers and a hidden dimension of 64 and a Trunk Network with four layers and a hidden dimension of 32. Across the considered values for these hyperparameters, this configuration gave the best result. For the Burgers equation, we used a Branch Network with a hidden dimension of 128 and a Trunk Network with a hidden dimension of 64. Still, the results were not competitive with NPR.To calculate the metrics for the heat equation, we compute a reference solution on a grid of size 500 2 using a finite difference scheme. Concretely, we discretize the domain [0, 1] × [0, 1] by using equidistant grid points x 1 , . . . , x n and t 1 , . . . , t m with m = n = 500.We then compute the absolute value of the difference d abs,ij between the reference solution and the model output on this grid and the errors asL 1 = 500 -2 xi,tj d abs,ij , L 2 = 500 -2 xi,tj d 2 abs,ij , L ∞ = max xi,tj d abs,ij .For the Burgers equation, we use the explicit formula from Equation 6and compute the errors in the same way as for the heat equation.The idea of PINNs (Raissi et al., 2019) is to design a loss function, which incorporates all the components of equation 1 and is zero if and only if the network solves the PDE. The typical approach is to parametrize the solution by an NN u θ with parameters θ and introduce three loss components:L PDE = [0,T ]×Ω ∥∂ t u θ (t, x) -N (u θ (t, x))∥ dxdt L IC = Ω ∥u θ (x, 0) -u 0 (x)∥ dx L BC = [0,T ]×∂Ω ∥B(u θ (t, x))∥ dxdt.The partial derivatives needed to evaluate ∂ t u θ and N u are computed using automatic differentiation and the norm ∥•∥ appearing in the integrals is typically the L2 norm. The total loss is given by L = λ PDE L PDE + λ IC L IC + λ BC L BC , where λ PDE , λ IC and λ BC are loss weights which can be chosen as hyperparameters or adapted dynamically (Wang et al., 2021a; 2022b) .The DeepONet (Lu et al., 2021a) architecture gives an approximation of any operator G implicitly, that means: Given a function u and a value y = (t, x), the NN returns an approximation of G(u)(y).While this approach has shown to be successful in many applications, it has some drawbacks. First, the mapping u → G(u) is a bit hidden in the architecture. It can be recovered by evaluating the Branch Network G branch at u and then considering the mapping(t, x) → G trunk (t, x) • G branch (u).This type of recovery can be useful if one wants to evaluate the function G branch (u) at many points without having to run the full network. This however illustrates another shortcoming of the Deep-ONet architecture: An unneccesarily large amount of the complexity needs to be allocated at the Trunk Network, which must produce rich enough representations of the inputs (t, x) to be able to approximate G(u)(t, x) for any u by an operation as simple as the dot product, without knowledge of u. This issue has been addressed in Wang et al. (2022a) by introducing additional encoders which mix the hidden representations of the Branch Network and Trunk Network at every step. However, introducing these encoders makes the recovery of the mapping u → G(u) impossible, as now for each pair (u, (t, x)) the whole architecture needs to be evaluated.Building upon the PINN framework, the Physics-informed DeepONet (Wang et al., 2021b) applies these principles in an operator learning context. It harnesses the power of automatic differentiation and extends the capabilities of PINNs to learn mappings between function spaces while respecting physical constraints. It is worth emphasizing that both approaches do not require any labeled training data (except for initial-and boundary conditions, which are needed to make the problem well-posed) and can thus be considered as self-supervised learning methods. This is particularly crucial in the context of learning solution operators, as obtaining training data involves applying numerical solvers. Recently, Hasani & Ward (2024) have proposed to generate training data for neural operators which avoids solving the PDE numerically.In the specific setting of learning the mapping from an initial condition to the solution at a given time, a model N θ is posed with the challenge of learning an identity mapping for t = 0, i.e. N θ (u 0 , 0, x) = u 0 (x) must hold for all x ∈ Ω. Lu et al. (2021b) proposed to parametrize the output in a way that the initial condition is always satisfied. Concretely, given any model N θ , a new model Nθ is defined as follows:Nθ (u 0 , t, x) = (1 -α(t))N θ (u 0 , t, x) + α(t)u 0 (x).Here, α : [0, T ] → [0, 1] is a function which is 1 for t = 0 and 0 for t = T , e.g α(t) = 1 -t T . 2 This was adopted by Brecht et al. (2023) in the context of operator learning. This way, we do not only mitigate the problem of learning the identity, but also eliminate the need for the L IC term in the loss function, thus reducing the number of hyperparameters and making the problem easier. In the case of Dirichlet boundary conditions, we can apply the same procedure to those by parametrizing the output by a function β : [0, T ] × ∂Ω → R d and defining the new model as:EQUATIONwhere u b : [0, T ] → ∂Ω describes the (non-homogenuous) Dirichlet boundary condition and β : ∂Ω → [0, 1] is a function which is 1 for x ∈ ∂Ω and between 0 and 1 for x ∈ Ω. This eliminates L BC from the loss function. We point out that neither a closed form of u 0 nor u b is required. If we only have samples of u 0 and u b at certain points, we can simply interpolate them (e.g. using linear interpolation or splines).While we found that hardcoding the boundary conditions increases the speed at which our models converge, it is not strictly necessary. Hardcoding the initial conditions however is crucial. This is interesting, as it hints at the following: Regressing the parameters of an NN to act as a PINN is a very difficult task, but clearly it will be even harder if the networks are supposed to learn complicated mappings. Thus, our study demonstrates that hardcoding initial conditions greatly simplifies the learning problem and should be considered in any application of PINNs.A.7 HYPERNETWORKS Ha et al. (2017) have introduced the approach of one NN learning the parameters of another NN and have coined it Hypernetworks. We refer to Chauhan et al. ( 2023) for a recent review on the topic.A multilayer perceptron (MLP) is a NN with an input layer, an output layer and one or more hidden layers. Each layer consists of a linear transformation followed by a non-linear activation function.The linear transformation is parametrized by a matrix W and a bias vector b. Concretely, for an MLP with n hidden layers and hidden dimension d, the output is defined as:MLP(x) = W n-1 h n-1 + b n-1 h i = σ(W i h i-1 + b i ) for i = 2, . . . , n -1 h 1 = σ(W 1 x + b 1 )where x is the input, h i is the output of the i-th layer, W i are the weight matrices and b i are the bias vectors. σ is a nonlinear function, often called the activation and is typically chosen as tanh or ReLU.If we choose to parametrize the Target Network by a multilayer perceptron (MLP) with n hidden hidden layers, each with a hidden dimension of d hidden , the number of parameters is given byEQUATIONSince this scales quadratically in the hidden dimension d hidden , it quickly becomes intractable as it blows up the final layer of our network. Therefore, motivated by the success of LoRA (Hu et al., 2022; Zimmer et al., 2023a) , we parametrize the weights of the hidden layers as the product of 2 matrices, i.e. W i = A i B i for i = 1, . . . , n hidden -1 with A i , B T i ∈ R dhidden×r , where r ≪ d hidden is the rank of the matrices and therefore the rank of the matrix W i is also at most r. This reduces the number of parameters toEQUATIONRemarkably, for some PDEs we achieve good results with a rank as low as four, showing that by this approach we can learn approximations of the solution operator efficiently.Hypernetwork and Target Network Architectures Our Hypernetwork H Φ is parametrized as a multilayer perceptron (MLP). The choice of the number of hidden layers, their dimensions, and the activation function (sin instead of the common tanh) is made to allow for richer implicit neural representations, as suggested by Sitzmann et al. (2020) . The Target network T θ also follows an MLP architecture with low-rank hidden weights, designed to efficiently capture the dynamics of the problem.Sampling and Discretization of Initial Conditions Initial conditions u 0 are sampled from a compact set K ⊂ X using a specified sampling procedure. The discretization of u 0 into input vectors or tensors for H Φ in this work simply means taking equidistant samples of u 0 on the domain Ω. The number of samples is a hyperparameter of the model ans was chosen to be 32 in this study.Periodic Update of Loss Weights The loss weights λ PDE , λ IC , λ BC are updated periodically to normalize the magnitudes of the gradient updates. Considering the case that all loss components are present (neither initial-nor boundary conditions are hardcoded), every 100 steps we compute the loss weights as follows:1. Sample a batch for each of the loss components.2. Compute ∂Φ ∂Li for i in {PDE, IC, BC} using automatic differentiation.3. Compute the total magnitude of the updates M = i ∂Φ ∂Li . 
set λ: 
 i = M ∂Φ ∂Li -1for i in {PDE, IC, BC}.We explain in more detail how minibatches are constructed for computing the loss components. For L PDE , we sample n pde initial conditions from K as well as n pde space-time coordinates (t, x) from [0, T ] × Ω. The initial conditions u 0 are discretized into vectors û, so a batch which is input to the hypernetwork is of the form ûi , (t i , x i ), where i is the index for the position in the minibatch. After passing this batch through H Φ , we obtain a batch of vectors vi , which we reshape into a batch of NNs T θ,i . We then evaluate (in parallel) the NNs T θ,i at the points (t i , x i ) and use automatic differentiation to compute the derivatives ∂ t T θ,i and N T θ,i . Finally, we compute the loss component L PDE and the gradient of the hypernetwork parameters Φ w.r.t. L PDE using automatic differentiation. For L BC , we sample n bc initial conditions from K as well as n bc space-time coordinates from ∂Ω × [0, T ]. Again, the initial conditions are discretized and passed through the hypernetwork, after which we evaluate the NNs at the sampled points and compute the loss component L BC and its gradient in a supervised manner. For L IC , the batch consists of n ic initial conditions sampled from K together with points x i sampled from Ω. We pass the discretized initial conditions through the hypernetwork and evaluate the NNs at (0, x i ). We then compute the loss component L IC and its gradient in a supervised manner. Finally, we construct the total loss L = λ PDE L PDE + λ IC L IC + λ BC L BC and compute the gradient of the Hypernetwork parameters Φ w.r.t. L using automatic differentiation. We then update the Hypernetwork parameters using the Adam Algorithm. We present the full training algorithm in Algorithm 1.To enhance the adaptability of our model for both challenging in-distribution initial conditions and out-of-distribution examples, we propose a systematic procedure. Our model, by design, outputs NN parameters corresponding to each given initial condition. This capability enables us to input any initial condition, within the constraints of our sensor resolution3 , to obtain a specific set of NN parameters. To increase expressivity, we then "unfold" this network by explicitly computing the low-rank product matrices A i B i for the hidden weights. Subsequently, the Branch Network is omitted, treating the resulting network as a conventional PINN. Remarkably, this adaptation process, which requires only a few hundred steps or a matter of seconds, significantly reduces the error for a wide range of initial conditions, including out-of-distribution as shown in Figure 3 .We present the results for both equations in Table1. For the heat equation, the performance in terms of L 1 and L 2 errors is comparable to the DeepONet, the L ∞ error however is much lower across all considered configurations. This is likely due to the enhanced expressivity of NPR, which enables it to precisely capture the solution in the whole space-time domain. For the Burgers equation, NPR outperforms the DeepONet in all metrics and configurations, even though we increased the sizes of the Branch Network and Trunk Network as compared to the heat equation, underlining the superiority of NPR in capturing nonlinear dynamics. By design, NPR requires more parameters than the DeepONet, as the the output dimension of the Hypernetwork is a compressed representation of the parameters of a full NN. Increasing the amount of DeepONet parameters to a comparable amount gave worse results than the ones reported here.4 CONCLUSIONWe have presented Neural Parameter Regression, merging Hypernetworks and PINNs for selfsupervised learning of PDE solution operators. Utilizing low-rank matrices for compact parameterization and ensuring initial condition compliance, our approach efficiently learns PDE solution operators and adapts swiftly to new conditions, promising advancements in the field.While our method shows promising results, it shares a fundamental shortcoming of the DeepONet, in the sense that it is only applicable to operators between functions spaces of relatively low di-C k (Ω, R) denotes the k-times continuously differentiable functions from Ω to R.The factor α(t) in front of the u0 term is not necessary to achieve the hardcoding, however we found that it works better.This refers to initial conditions that can be meaningfully represented given the chosen number of sensors. 
 
 Appendix:  


Decision: This author presents a novel method for operator learning. Despite there are some concerns from reviewer YoR2, regarding the experiment results and whether this method can handle local features, those concerns can be address through rewriting. I strongly recommend the authors to address both reviewers' concerns in the camera-ready version.

Reviews: review: title: title: Interesting idea with good experimental evaluation.
main_review: main_review: The authors present a novel approach to operator learning with the application of learing solution operators for PDEs in mind. The proposed method applies the concept of Hypernetworks to the realm of Scientific Machine Learning and operator learning in particular.

The idea is simple and presented in a straight-forwards fashion. The idea is tested on two different PDE problems (Burgers equation and Heat equation) and hyperparameters are disclosed to make the setup reproducible. My main criticism is the limited evaluation of the method - as these two problems are fairly simple toy problems. As the self-supervised nature of physics-informed losses do not require extensive generation of training datasets, i would have welcomed more experiments, particularly in 2 spatial dimensions. Moreover, only DeepONets are considered as a baseline here. I understand that this is due to the problem setting. However, I would have welcomed a comparison to  other baselines such as Physics Informed Neural Operators (PINO).

Given that this is a workshop, I propose to accept the submission and encourage the authors to address my concern above, should they consider a full submission somewhere else.
title: title: Introduces Neural Parameter Regression (a model with two networks) for operator learning
main_review: main_review: The paper introduces Neural Parameter Regression or NPR for operator learning. The idea of the paper is to encode the imposed initial condition into a latent space which is then mapped to the solution space via a secondary neural network. The latter NN additionally takes (x, t) as inputs and avoids trivial solutions by reparametrizing the output such that the model learns the deviation from the initial conditions. 

I believe the paper has a novel idea but it is not analyzed with sufficient depth. I can see the following main issues with the current state of their approach: (1) in high spatial dimensions, how would you handle the large encoding space (i.e., the method does not seem to scale well to cases where the initial condition is specified in 3D), (2) any version of the model considered in Table 1 has much more parameters than DeepONet so I am not confused about its higher performance (i.e., a fair comparison would assign roughly the same number of parameters to each model), (3) is the method able to handle local features (e.g., high frequency or large gradients)? The examples used are very simple. (4) the method is discretization dependents (due to discretization of u0) so perhaps you can use PCA to remove this limitation.