Learning nonlinear level sets for dimensionality reduction in function approximation
Guannan Zhang Jiaxin Zhang Jacob Hinkle



 
 



Wide band spectroscopic response of monocrystallines to low dose neutron and gamma radiation
    [
    December 30, 2023
============================================================================================






We developed a Nonlinear Level-set Learning (NLL) method for dimensionality reduction in high-dimensional function approximation with small data. This work is motivated by a variety of design tasks in real-world engineering applications, where practitioners would replace their computationally intensive physical models (e.g., high-resolution fluid simulators) with fast-to-evaluate predictive machine learning models, so as to accelerate the engineering design processes. There are two major challenges in constructing such predictive models: (a) high-dimensional inputs (e.g., many independent design parameters) and (b) small training data, generated by running extremely time-consuming simulations. Thus, reducing the input dimension is critical to alleviate the over-fitting issue caused by data insufficiency. Existing methods, including sliced inverse regression and active subspace approaches, reduce the input dimension by learning a linear coordinate transformation; our main contribution is to extend the transformation approach to a nonlinear regime. Specifically, we exploit reversible networks (RevNets) to learn nonlinear level sets of a high-dimensional function and parameterize its level sets in low-dimensional spaces. A new loss function was designed to utilize samples of the target functions' gradient to encourage the transformed function to be sensitive to only a few transformed coordinates.
The NLL approach is demonstrated by applying it to three 2D functions and two 20D functions for showing the improved approximation accuracy with the use of nonlinear transformation, as well as to an 8D composite material design problem for optimizing the buckling-resistance performance of composite shells of rocket inter-stages.






§ INTRODUCTION



High-dimensional function approximation arises in a variety of engineering applications where scientists or engineers rely on accurate and fast-to-evaluate approximators to replace complex and time-consuming physical models (e.g. multiscale fluid models), so as to accelerate scientific discovery or engineering design/manufacture. In most of those applications, training and validation data need to be generated by running expensive simulations,

that the amount of training data is often limited due to high cost of data generation (see <ref> for an example). Thus, this effort is motivated by the challenge imposed by high dimensionality and small data in the context of function approximation. 




One way to overcome the challenge is to develop dimensionality reduction methods that can build a transformation of the input space to increase the anisotropy of the input-output map. In this work, we assume that the function has a scalar output and the input consists of high-dimensional independent variables, such that there is no intrinsically low-dimensional structure of the input manifold. In this case, instead of analyzing the input or output manifold separately, we will learn low-dimensional structures of the target function's level sets to reduce the input dimension. 

Several methods have been developed for this purpose, including sliced inverse regression and active subspace methods. A literature review of those methods is given in <ref>.

Despite many successful applications of those methods, 
their main drawback is that they use linear transformations to capture low-dimensional structures of level sets. When the geometry of the level sets is highly nonlinear, e.g., f( x) = sin( x_2^2), the number of active input dimensions cannot be reduced by linear transformations. 




In this effort, we exploited reversible residual neural networks (RevNets) <cit.> to learn the
target functions' level sets and build nonlinear coordinate transformations to reduce the number of active input dimensions of the function.
Reversible architectures have been developed in the literature <cit.> with a purpose of reducing memory usage in backward propagation, while we intend to exploit the reversibility to build bijective nonlinear transformations. Since the RevNet is used for a different purpose, we designed a new loss function for training the RevNets, such that a well-trained RevNet can capture the nonlinear geometry of the level sets. The key idea is to utilize samples of the function's gradient to promote the objective that most of the transformed coordinates move on the tangent planes of the target function, i.e., the transformed function is invariant with respect to those coordinates. In addition, we constrain the determinant of the Jacobian matrix of the transformation in order to enforce invertibility. 

The main contributions of this effort can be summarized as follows:



   
  * Development of RevNet-based coordinate transformation model for capturing the geometry of level sets, which extends function dimensionality reduction to the nonlinear regime. 
    
    
  * Design of a new loss function that exploits gradient of the target function to successfully train the proposed RevNet-based nonlinear transformation. 
    
    
  * Demonstration of the performance of the proposed NLL method on a high-dimensional real-world composite material design problem for rocket inter-stage manufacture.

 






















§ PROBLEM FORMULATION


We are interested in approximating a d-dimensional multivariate function of the form

    y = f( x),    x ∈Ω⊂ℝ^d,

where Ω is a bounded domain in ℝ^d, the input x := (x_1, x_2, …, x_d)^⊤ is a d-dimensional vector, and the output y is a scalar value. Ω is equipped with a probability density function ρ: ℝ^d ↦ℝ^+, i.e.,

    0 < ρ( x) < ∞ ,   x ∈Ω  and  ρ( x) = 0,   x ∉Ω,

and all the components of x are assumed to be independent. The target function f is assumed to be first-order continuously differentiable, i.e., f ∈ C^1(Ω), and square-integrable with respect to the probability measure ρ, i.e., ∫_Ω f^2( x) ρ( x) d x < ∞. 

In many engineering applications, e.g., the composite shell design problem in <ref>, f usually represents the input-output relationship of computationally expensive simulators. In order to accelerate a discovery/design process, practitioners seek to build an approximation of f, denoted by f̃, such that the error f-f̃ is smaller than a prescribed threshold ε > 0, i.e.,
f( x)-f̃( x)_L^2_ρ(Ω) < ε,
where ·_L^2_ρ is the L^2 norm under the probability measure ρ. 

As discussed in <ref>, the main challenge results from the concurrence of having high-dimensional input (i.e., large d) and small training data, which means the amount of training data is insufficient to overcome the curse of dimensionality. In this scenario, naive applications of existing approximation methods, e.g., sparse polynomials, kernel methods, neural networks (NN), etc., may lead to severe over-fitting. Therefore, our goal is to reduce the input dimension d by transforming the original input vector x to a lower-dimensional vector z, such that the transformed function can be accurately approximated with small data. 





 §.§ Related work







Manifold learning for dimensionality reduction. Manifold learning, including linear and nonlinear approaches <cit.>, focuses on reducing data dimension via learning intrinsically low-dimensional structures in the data. Nevertheless, since we assume the input vector x in Eq. (<ref>) consists of independent components and the output f is a scalar, no low-dimensional structure can be identified by separately analyzing the input and the output data. Thus, the standard manifold learning approaches are not applicable to the function dimensionality reduction problem under consideration.




Sliced inverse regression (SIR). SIR is a statistical dimensionality reduction approach for the problem under consideration. In SIR, the input dimension is reduced by constructing/learning a linear coordinate transformation z = 𝐀 x, with the expectation that the output of the transformed function y = h( z ) = h(𝐀 x) is only sensitive to a very small number of the new coordinates of z. The original version of SIR was developed in <cit.> and then improved extensively by <cit.>. To relax the elliptic assumption (e.g., Gaussian) of the data, kernel dimension reduction (KDR) was introduced in <cit.>. 

Several recent work, including manifold learning with KDR <cit.> and localized SIR <cit.>, were developed for classification problem. In <ref>, the SIR will be used to produce baseline results to compare with the performance of our nonlinear method.







Active subspace (AS). The AS method <cit.> shares the same motivation as SIR, i.e., reducing the input dimension of f( x) by defining a linear transformation z = 𝐀 x. The main difference between AS and SIR is the way to construct the matrix 𝐀.
The AS method does not need the elliptic assumption needed for SIR but requires (approximate) gradient samples of f( x) to build 𝐀. For both SIR and AS, when the level sets of f are nonlinear, e.g., f( x) = sin( x^2_2), the dimension cannot be effectively reduced using any linear transformation. The AS method will be used as another baseline to compare with our method in <ref>.




Reversible neural networks. We exploited the RevNets proposed in <cit.> to define our nonlinear transformation for dimensionality reduction. Those RevNets describe bijective continuous dynamics while regular residual networks result in crossing and collapsing paths which correspond to non-bijective continuous dynamics <cit.>. Recently, 
RevNets have been shown to produce competitive performance on discriminative tasks <cit.> and generative tasks <cit.>. In particular,
the nonlinear independent component estimation (NICE) <cit.> used RevNets to build nonlinear coordinate transformations to factorize high-dimensional density functions into products of independent 1D distributions. The main difference between NICE and our approach is that NICE seeks convergence in distribution (weak convergence) with the purpose of building an easy-to-sample distribution, and our approach seeks strong convergence as indicated by the norm ·_L_ρ^2 with the purpose of building an accurate pointwise approximation to the target function in a lower-dimensional input space.













§ PROPOSED METHOD: NONLINEAR LEVEL SETS LEARNING (NLL)

The goal of dimensionality reduction is 
to construct a bijective nonlinear transformation, denoted by

    z =  g( x) ∈ℝ^d   and   x =  g^-1( z),

where z = (z_1, …, z_d)^⊤, such that the composite function
y = f ∘ g^-1 ( z)
has a very small number of active input components. In other words, even though z ∈ℝ^d is still defined in ℝ^d, the components of z
can be split into two groups, i.e.,
z = ( z_ act,  z_ inact) 
with dim( z_ act) much smaller than d, such that f∘ g^-1 is only sensitive to the perturbation of z_ act. To this end, our method was inspired by the following observation:

Observation: For a fixed pair ( x,  z) satisfying z =  g( x), if x =  g^-1( z), as a particle in Ω, moves along a tangent direction of the level set passing through f( x) under a perturbation of z_i (the i-th component of z), then the output of f∘ g^-1 ( z) does NOT change with z_i in the neighbourhood of z.


Based on such observation, we intend to build and train a nonlinear transformation g with the objective that having a prescribed number of inactive components of z satisfy the above statement, and those inactive components will form z_ inact.  






Training data: 
We need two types of data for training g, i.e., samples of the function values and its gradients, denoted by

    Ξ_ train := {( x^(s), f( x^(s)), ∇ f( x^(s)) ):   s = 1, …, S },

where { x^(s):s=1, …, S} are drawn from ρ( x), and ∇ f( x^(s)) denotes the gradient of f at x^(s). The gradient samples describe the tangent direction of the target function's level sets, i.e., the gradient direction is in perpendicular to all tangent directions. The requirement of gradient samples may limit the applicability of our approach to real-world applications in which gradient is not available. A detailed discussion on how to mitigate such disadvantage is given in <ref>.






 §.§ The level sets learning model: RevNets

The first step is to define a model for the nonlinear transformation g in Eq. (<ref>). In this effort, we utilize the nonlinear RevNet model proposed in <cit.>, defined by













    {    u_n+1 =  u_n + h 𝐊_n,1^⊤ σ(𝐊_n,1 v_n +  b_n,1),
        v_n+1 =  v_n  - h 𝐊_n,2^⊤ σ(𝐊_n,2 u_n+1 +  b_n,2),
    .

for n = 0, 1, …, N-1, where u_n and v_n are partitions of the states, h is the “time step” scalar, 𝐊_n,1, 𝐊_n,2 are weight matrices, b_n,1, b_n,2 are biases, and σ is the activation function.
Since u_n, v_n can be explicitly calculated given u_n+1,  v_n+1, the RevNet in Eq. (<ref>) is reversible by definition. Even though our approach can incorporate any reversible architecture, we chose the model in Eq. (<ref>) because it has been shown in <cit.> that this architecture has better nonlinear representability than other types of RevNets. 

To define g:  x ↦ z, we split the components of x evenly into 
u_0 and v_0, and split the components of z accordingly into u_N and v_N, i.e.

    x    := 
    [ u_0; v_0 ]  where    u_0 := (x_1, …, x_⌈ d/2 ⌉)^⊤,  v_0 := (x_⌈ d/2 ⌉+1, …, x_d)^⊤, 
      z    := 
    [ u_N; v_N ]  where    u_N := (z_1, …, z_⌈ d/2 ⌉)^⊤,  v_N := (z_⌈ d/2 ⌉+1, …, z_d)^⊤,

such that the nonlinear transformation g is defined by the map ( u_0,  v_0) ↦ ( u_N,  v_N) from the input states of the N-layer RevNets in Eq. (<ref>) to its output states, i.e.,


    x = [ u_0; v_0 ]
    
    [ u_N; v_N ] =  z.

It was shown in <cit.> that the RevNet in Eq. (<ref>) is guaranteed to be stable, so that we can use deep architectures to build a highly nonlinear transformation to capture the geometry of the level sets of f.




 §.§ The loss function

The main novelty of this work is the design of the loss function for training the RevNet in Eq. (<ref>). The new loss function includes two components. The first component is to 
inspired by our observation given at the beginning of <ref>. Specifically, 





guided by such observation, we write out the Jacobian matrix of the inverse transformation g^-1:  z ↦ x as



    𝐉_g^-1( z) = [ J_1( z),  J_2( z), …,  J_d(z)]   with   J_i( z) := (∂ x_1∂ z_i( z), …, ∂ x_d∂ z_i( z))^⊤


where the i-th column J_i describes the direction in which the particle x moves when perturbing z_i. As such, we can use J_i to mathematically rewrite our observation as: the output of f( x) does not change with a perturbation of z_i in the neighborhood of z, if


    J_i( z) ∇ f( x) ⟺ ⟨ J_i( z),  ∇ f( x) ⟩= 0,


where ⟨·, ·⟩ denotes the inner product of two vectors. 
The relation in Eq. (<ref>) is illustrated in Figure <ref>.








Therefore, the first component of the loss function, denoted by L_1, is defined by

    L_1 := ∑_s=1^S ∑_i=1^d [ω_i ⟨ J_i( z^(s)) J_i( z^(s))_2,  ∇ f( x^(s)) ⟩]^2,

where ω_1, ω_2, …, ω_d are user-defined anisotropy weights determining how strict the condition in Eq. (<ref>) is enforced on each dimension. A extreme case could be ω := (0, 1,1,…, 1), which means the objective is to train the transformation g such that the intrinsic dimension of f∘ g^-1( z) is one when L_1 = 0. Another extreme case is ω = (0, …, 0), which leads to L_1 = 0 and no dimensionality reduction will be performed. In practice, such weights ω give us the flexibility to balance between training cost and reduction effect. 

It should be noted that we only normalize J_i in Eq (<ref>), but not ∇ f, such that L_1 will not penalize too much in the regions where ∇ f is very small. 
In particular, L_1 = 0 if f is a constant function. 

The second component of the loss function is designed to guarantee that the nonlinear transformation
r0.3


  


    < g r a p h i c s >

  

  

 Illustration of the observation for defining the loss L_1 in Eq. (<ref>), i.e., f( x) is insensitive to perturbation of z_i in the neighborhood of z if J_i( z) ∇ f( x), where J_i is defined in Eq. (<ref>).



 g is non-singular. It is observed in Eq. (<ref>) that L_1 only affects the Jacobian columns J_i with ω_i ≠ 0, but has no control of the columns J_i with ω_i = 0. To avoid the situation that the transformation g becomes singular during training, we define the second loss component L_2 as a quadratic penalty on the Jacobian determinant, i.e. ,




    L_2 := (  det(𝐉_g^-1) - 1)^2,

which will push the transformation to be non-singular and volume preserving. Note that L_2 can be viewed as a regularization term.   














In summary, the final loss function is defined by


    L := L_1 + λ L_2,

where λ is a user-specified constant to balance the two terms.



 §.§ Implementation

The RevNet in Eq. (<ref>) with the new loss function in Eq. (<ref>) was implemented in PyTorch 1.1 and tested on a 2014 iMac Desktop with a 4 GHz Intel Core i7 CPU and 32 GB DDR3 memory. To make use of the automatic differentiation in PyTorch, we implemented a customized loss function in Pytorch, where the entries of the Jacobian matrix 𝐉_g^-1 were computed using finite difference schemes, and the Jacobian det(𝐉_g^-1) was approximately calculated using the PyTorch version of singular value decomposition. Since this effort focuses on proof of concept of the proposed methodology, the current implementation is not optimized in terms of computational efficiency.




§ NUMERICAL EXPERIMENTS

We evaluated our method using three 2D functions in <ref> for visualizing the nonlinear capability, two 20D functions in <ref> for comparing our method with brute-force neural networks, SIR and AS methods, as well as a composite material design problem in <ref> for demonstrating the potential impact of our method on real-world engineering problems. To generate baseline results, we used existing SIR and AS codes available at <https://github.com/paulcon/active_subspaces> and <https://github.com/joshloyal/sliced>, respectively. Source code for the proposed NLL method is available in the supplemental material.  










 §.§ Tests on two-dimensional functions

Here we applied our method to the following three 2-dimensional functions: 


    f_1( x)    = 1/2sin(2π( x_1 + x_2)) +1   for  x ∈Ω = [0,1]×[0,1] 
    
    
    f_2( x)    = exp(- (x_1-0.5)^2 - x_2^2)   for 
    x ∈Ω = [0,1]×[0,1] 
    
    
    f_3( x)    = x_1^3 + x_2^3 + 0.2 x_1 + 0.6 x_2   for  x ∈Ω = [-1,1]×[-1,1].

We used the same RevNet architecture for the three functions. Specifically, u and v in Eq. (<ref>) were 1D variables (as the total dimension is 2); the number of layers was N = 10, i.e., 10 blocks of the form in Eq. (<ref>) were connected; 𝐊_n,1, 𝐊_n,2 were 2× 1 matrices; b_n,1,  b_n,2 are 2D vectors; the activation function was tanh; the time step h was set to 0.25; stochastic gradient descent method was used to train the RevNet with the learning rate being 0.01; no regularization was applied to the network parameters; the weights in Eq. (<ref>) was set to ω = (0,1); λ = 1 in the loss function in Eq. (<ref>); the training set included 121 uniformly distributed samples in Ω, and the validation set included 2000 uniformly distributed samples in Ω. We compared our method with either SIR or AS for each of the three functions.

The results for f_1, f_2, f_3 are shown in Figure <ref>. For f_1, it is known that the optimal transformation is a 45 degree rotation of the original coordinate system.
The first row in Figure <ref> shows that the trained RevNet can approximately recover the 45 degree rotation, which demonstrates that the NLL method can also recover linear transformation. The level sets of f_2 and f_3 are nonlinear, and the NLL method successfully captured such nonlinearity. In comparison, the performance of AS and SIR is worse than the NLL method because they can only perform linear transformation.








 §.§ Tests on 20-dimensional functions








Here we applied the new method to the following two 20-dimensional functions: 




    f_4( x) = sin(x_1^2 + x_2^2 + ⋯ + x_20^2)   and    f_5( x) = ∏_i=1^20(1.2^-2 + x_i^2 )^-1





for x ∈Ω = [0,1]^20.We used one RevNet architecture for the two functions. Specifically, u and v in Eq. (<ref>) were 10D variables, respectively; the number of layers is N = 30, i.e., 30 blocks of the form in Eq. (<ref>) were connected; 𝐊_n,1, 𝐊_n,2 were 20× 10 matrices; b_n,1,  b_n,2 were 10-dimensional vectors; the activation function was tanh; the time step h was set to 0.25; stochastic gradient descent method was used to train the RevNet with the learning rate being 0.05; λ = 1 for the loss function in Eq. (<ref>); the training set includes 500 uniformly distributed samples in Ω.




























The effectiveness of the NLL method is shown as relative sensitivity indicators in Figure <ref>(a) for f_4 and Figure <ref>(b) for f_5. 
The sensitivity of each transformed variable z_i is described by the normalized sample mean of the absolute values of the corresponding partial derivative. The definition of w in Eq. (<ref>) provides the target anisotropy of the transformed function. For f_4, we set ω_1 = 0, ω_i = 1 for i = 2, …, 20; for f_5 we set ω_1 = ω_2 = 0, ω_i = 1 for i = 3, …, 20. As expected, the NLL method successfully reduced the sensitivities of the inactive dimensions to two orders of magnitude smaller than the active dimensions. In comparison, the SIR and AS methods can only reduce their sensitivities by one order of magnitude using optimal linear transformations.

Next, we show how the NLL method improves the accuracy of the approximation of the transformed function f∘ g^-1. We used two fully connected NNs to approximate the transformed functions, i.e., one has 2 hidden layers with 20+20 neurons, and the other has a single hidden layer with 10 neurons. The implementation of both networks was based on the neural network toolbox in Matlab 2017a. We used various sizes of training data: 100, 500, 10,000, and we used another 10,000 samples as validation data. All the samples are drawn uniformly in Ω.
The approximation error was computed as the relative root mean square error (RMSE) using the validation data. For comparison, we used the same data to run brute-force neural networks without any transformation, AS and SIR methods. 

The results for f_4 and f_5 are shown in Table <ref> and <ref>, respectively. For the 20+20 network, when the training data is too small (e.g., 100 samples), all the methods have the over-fitting issue; when the training data is very large (e.g., 10,000 samples), all the methods can achieve good accuracy
[ 10% or smaller RMSE is considered as satisfactory accuracy in many engineering applications]. Our method shows significant advantages over AS and SIR methods, when having relatively small training data, e.g., 500 training data, which is a common scenario in scientific and engineering applications. For the single hidden layer network with 10 neurons, we can see that the brute-force NN, AS and SIR cannot achieve good accuracy with 10,000 training data (no over-fitting), which means the network does not have sufficient expressive power to approximate the original function and the transformed functions using AS or SIR. In comparison, the NLL method still performs well as shown in Table <ref>(Right) and <ref>(Right). This means the dimensionality reduction has significantly simplified the target functions' structure, such that the transformed functions can be accurately approximated with smaller architectures to reduce the possibility of over-fitting.

























  
  
  
  
    
    
      
      
     
        
    
    
    
    
    
    
    
    
     
      
      
     
        
    
    
    
    
    
    
    
  
  
  
  
   
    
    
      
      
    
        
    
    
    
    
    
    
    
    
      
      
    
        
    
    
    
    
    
    
    
  
















































  
  
  
  
    
    
      
      
     
        
    
    
    
    
    
    
    
     
      
      
     
        
    
    
    
    
    
    
    
  
  
  
  
   
    
    
      
      
    
        
    
    
    
    
    
    
    
    
      
      
    
        
    
    
    
    
    
    
    
  






















































































 §.§ Design of composite shell for rocket inter-stages

Finally, we demonstrate the NLL method on a real-world composite material design problem. 
With high specific stiffness and strength, composite materials are increasingly being used for launch-vehicle structures. A series of large-scale composite tests for shell buckling 
knockdown factor conducted by NASA (see Figure <ref>(a)) aimed to develop and validate new analysis-based design guidelines for safer and lighter space structure <cit.>. Since the experimental cost is extremely high, numerical simulation, e.g., finite element method (FEM), is often employed to predict the shell buckling knockdown factor given a multi-layer ply stack design<cit.>, as illustrated in Figure <ref>(c). The goal of this work is to implement an accurate approximation of this high-dimensional regression problem where the inputs are ply angles for 8 layers and the output is the knockdown factor which needs high precision for space structure design. However, the high fidelity FEM simulation is so time consuming that one analysis takes 10 hours and consequently it is impractical to collect a large data set for approximating the knockdown factor. 


To demonstrate the applicability of our method to this problem, we used a simplified FEM model
that runs relatively faster but preserves all the physical properties as the high-fidelity FEM model. As shown in Figure <ref>(b), a ply angle ranging from 0^∘ to 22.5^∘ that is assigned for each of 8 layers are considered in this example, i.e., the input domain is Ω = [0^∘,22.5^∘]^8. The RevNet has N = 10 layers; 𝐊_n,1, 𝐊_n,2 were 8× 4 matrices; b_n,1,  b_n,2 were 4-dimensional vectors; the activation function was tanh; the time step h was set to 0.1; stochastic gradient descent method was used with the learning rate being 0.05; λ = 1 for the loss function in Eq. (<ref>).




Like previous examples, we show the comparison of relative sensitivities in Table <ref>, where we allowed 3 active dimensions in the loss L_1, i.e. ω_1 = ω_2 = ω_3 = 0, and ω_i = 1 for i = 4,…, 8. 
r6.5cm
 

  Relative RMSE for approximating the composite material design model.
  

  
        2c100 data    
      2c500 data     
      2c10,000 data 

     (r)2-3 (r)4-5 (r)6-7
            Valid      Train     Valid     Train     Valid     Train

    NN             65.74%      0.01%     67.57%     24.77%     3.74%      3.52% 

    [gray]0.8                                                  
    NLL            63.18%      0.02%     11.96%     5.13%      2.51%      2.17% 

    AS             58.89%      0.13%     47.27%     19.11%     3.05%      2.91% 

    SIR            65.34%      0.21%     54.99%     22.52%     3.32%      3.21% 

    
  

  

As expected, the NLL method successfully reduced the input dimension by reducing the sensitivities of Dim 4-8 to two orders of magnitude smaller than the most active dimension, which outperforms the AS and SIR method. In Table <ref>, we show the RMSE of approximating the transformed function using a neural network with a single hidden layer having 20 neurons. The other settings are the same as the examples in <ref>. As expected, the NLL approach outperforms the AS and SIR in the small data regime, i.e., 500 training data. In Figure <ref>, we show the decay of the loss function for different choices of the anisotropy weights ω in L_1, we can see that the more inactive/insensitive dimensions (more non-zero ω_i), the slower the loss function decay.








  
  
    
    
      
      
     
        
    
    
    
    
    
    
    
     
      
      
     
        
    
    
    
    
    
    
    
  
  
  
  
  
  





  
  



























































































































































































§ CONCLUDING REMARKS

We developed RevNet-based level sets learning method for dimensionality reduction in high-dimensional function approximation. With a custom-designed loss function, the RevNet-based nonlinear transformation can effectively learn the nonlinearity of the target function's level sets, so that the input dimension can be significantly reduced. 

Limitations. Despite the successful applications of the NLL method shown in <ref>, we realize that there are several limitations with the current NLL algorithm, including (a) The need for gradient samples. Many engineering models do not provide gradient as an output. To use the current algorithm, we need to compute the gradients by finite difference or other perturbation methods, which will increase the computational cost. (b) Non-uniqueness. Unlike the AS and SIR method, the nonlinear transformation produced by the NLL method is not unique, which poses a challenge in the design of the RevNet architectures. (c) High cost of computing Jacobians. The main cost in the backward propagation lies in the computation of the Jacobian matrices and its determinant, which deteriorates the training efficiency as we increase the depth of the RevNet.

Future work.
There are several research directions we will pursue in the future. The first is to develop a gradient estimation approach that can approximately compute gradients needed by our approach. Specifically, we will exploit the contour regression method <cit.> and the manifold tangent learning approach <cit.>, both of which have the potential to estimate gradients by using function samples. The second is to improve the computational efficiency of the training algorithm. Since our loss function is more complicated than standard loss functions, it will require extra effort to improve the efficiency of backward propagation. 








99

behrmann2018invertible
Jens Behrmann, David Duvenaud, and Jörn-Henrik Jacobsen.
Invertible residual networks.
arXiv preprint arXiv:1811.00995, 2018.

belkin2006manifold
Mikhail Belkin, Partha Niyogi, and Vikas Sindhwani.
Manifold regularization: A geometric framework for learning from
  labeled and unlabeled examples.
Journal of machine learning research, 7(Nov):2399–2434, 2006.

NIPS2004_2647
Yoshua Bengio and Martin Monperrus.
Non-local manifold tangent learning.
In L. K. Saul, Y. Weiss, and L. Bottou, editors, Advances in
  Neural Information Processing Systems 17, pages 129–136. MIT Press, 2005.

castro2014geometric
Saullo GP Castro, Rolf Zimmermann, Mariano A Arbelo, Regina Khakimova, Mark W
  Hilburger, and Richard Degenhardt.
Geometric imperfections and lower-bound methods used to calculate
  knock-down factors for axially compressed composite cylindrical shells.
Thin-Walled Structures, 74:118–132, 2014.

AAAI1816517
Bo Chang, Lili Meng, Eldad Haber, Lars Ruthotto, David Begert, and Elliot
  Holtham.
Reversible architectures for arbitrarily deep residual neural
  networks.
In AAAI Conference on Artificial Intelligence, 2018.

constantine2015active
Paul G Constantine.
Active subspaces: Emerging ideas for dimension reduction in
  parameter studies, volume 2.
SIAM, 2015.

constantine2014active
Paul G Constantine, Eric Dow, and Qiqi Wang.
Active subspace methods in theory and practice: applications to
  kriging surfaces.
SIAM Journal on Scientific Computing, 36(4):A1500–A1524, 2014.

cook2005sufficient
R Dennis Cook and Liqiang Ni.
Sufficient dimension reduction via inverse regression: A minimum
  discrepancy approach.
Journal of the American Statistical Association,
  100(470):410–428, 2005.

cook1991sliced
R Dennis Cook and Sanford Weisberg.
Sliced inverse regression for dimension reduction: Comment.
Journal of the American Statistical Association,
  86(414):328–332, 1991.

cook2001theory
R Dennis Cook and Xiangrong Yin.
Theory & methods: special invited paper: dimension reduction and
  visualization in discriminant analysis (with discussion).
Australian & New Zealand Journal of Statistics,
  43(2):147–199, 2001.

dinh2014nice
Laurent Dinh, David Krueger, and Yoshua Bengio.
Nice: Non-linear independent components estimation.
arXiv preprint arXiv:1410.8516, 2014.

dinh2016density
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.
Density estimation using real nvp.
arXiv preprint arXiv:1605.08803, 2016.

donoho2003hessian
David L Donoho and Carrie Grimes.
Hessian eigenmaps: Locally linear embedding techniques for
  high-dimensional data.
Proceedings of the National Academy of Sciences,
  100(10):5591–5596, 2003.

fukumizu2004dimensionality
Kenji Fukumizu, Francis R Bach, and Michael I Jordan.
Dimensionality reduction for supervised learning with reproducing
  kernel hilbert spaces.
Journal of Machine Learning Research, 5(Jan):73–99, 2004.

fukumizu2009kernel
Kenji Fukumizu, Francis R Bach, Michael I Jordan, et al.
Kernel dimension reduction in regression.
The Annals of Statistics, 37(4):1871–1905, 2009.

gomez2017reversible
Aidan N Gomez, Mengye Ren, Raquel Urtasun, and Roger B Grosse.
The reversible residual network: Backpropagation without storing
  activations.
In Advances in neural information processing systems, pages
  2214–2224, 2017.

2018InvPr..34a4004H
Eldad Haber and Lars Ruthotto.
Stable architectures for deep neural networks.
Inverse Problems, 34:014004, Jan 2018.

Hauser:2017ta
Michael Hauser and Asok Ray.
Principles of Riemannian Geometry in Neural Networks.
NIPS, 2017.

hilburger2015shell
Mark W Hilburger.
Shell buckling knockdown factor project overview and status.
2015.

jacobsen2018revnet
Jörn-Henrik Jacobsen, Arnold Smeulders, and Edouard Oyallon.
i-revnet: Deep invertible networks.
arXiv preprint arXiv:1802.07088, 2018.

kingma2018glow
Durk P Kingma and Prafulla Dhariwal.
Glow: Generative flow with invertible 1x1 convolutions.
In Advances in Neural Information Processing Systems, pages
  10215–10224, 2018.

li2005
Bing Li, Hongyuan Zha, and Francesca Chiaromonte.
Contour regression: A general approach to dimension reduction.
Ann. Statist., 33(4):1580–1616, 08 2005.

li1991sliced
Ker-Chau Li.
Sliced inverse regression for dimension reduction.
Journal of the American Statistical Association,
  86(414):316–327, 1991.

li1992principal
Ker-Chau Li.
On principal hessian directions for data visualization and dimension
  reduction: Another application of stein's lemma.
Journal of the American Statistical Association,
  87(420):1025–1039, 1992.

li2007sparse
Lexin Li.
Sparse sufficient dimension reduction.
Biometrika, 94(3):603–613, 2007.

park2015bayesian
Mijung Park, Wittawat Jitkrittum, Ahmad Qamar, Zoltán Szabó, Lars
  Buesing, and Maneesh Sahani.
Bayesian manifold learning: the locally linear latent variable model
  (ll-lvm).
In Advances in neural information processing systems, pages
  154–162, 2015.

roweis2000nonlinear
Sam T Roweis and Lawrence K Saul.
Nonlinear dimensionality reduction by locally linear embedding.
science, 290(5500):2323–2326, 2000.

tenenbaum2000global
Joshua B Tenenbaum, Vin De Silva, and John C Langford.
A global geometric framework for nonlinear dimensionality reduction.
science, 290(5500):2319–2323, 2000.

wang2005adaptive
Jing Wang, Zhenyue Zhang, and Hongyuan Zha.
Adaptive manifold learning.
In Advances in neural information processing systems, pages
  1473–1480, 2005.

wu2009localized
Qiang Wu, Sayan Mukherjee, and Feng Liang.
Localized sliced inverse regression.
In Advances in neural information processing systems, pages
  1785–1792, 2009.

yeh2008nonlinear
Yi-Ren Yeh, Su-Yun Huang, and Yuh-Jye Lee.
Nonlinear dimension reduction with kernel sliced inverse regression.
IEEE transactions on Knowledge and Data Engineering,
  21(11):1590–1603, 2008.









