A VARIABLE P ROJECTION M ETHOD FOR C OMPUTATIONAL PDE SWITH A RTIFICIAL N EURAL N ETWORKS

Suchuan Dong
Center for Computational & Applied Mathematics
Department of Mathematics
Purdue University
West Lafayette, IN 47907
sdong@purdue.edu

A BSTRACTWe present a method for solving linear and nonlinear partial differential equations (PDE) basedon the variable projection framework and artificial neural networks. For linear PDEs, enforcing
the boundary/initial value problem on the collocation points gives rise to a separable nonlinear
least squares problem about the network coefficients. We reformulate this problem by the variable
projection approach to eliminate the linear output-layer coefficients, leading to a reduced problem
about the hidden-layer coefficients only. The reduced problem is computed first by the nonlinear leastsquares method to determine the hidden-layer coefficients, and then the output-layer coefficients are
determined by the linear least squares method. For nonlinear PDEs, enforcing the boundary/initialvalue problem on the collocation points gives rise to a nonlinear least squares problem that is notseparable. To enable the variable projection approach for nonlinear PDEs, we first linearize theproblem with a Newton iteration. The linearized system is solved by the variable projection framework
together with artificial neural networks. The current method exhibits a spectral-like accuracy, with its
errors decreasing exponentially with respect to the number of collocation points or the output-layer
coefficients. We compare this method with PINN, and show that the current method is much moreaccurate and computationally efficient.
Keywords Variable projection · Neural Network · Nolinear least squares · Spectral accuracy

1

Introduction

This work concerns the approximation of partial differential equations (PDE) with artificial neural networks (ANN), andwe exploit the the variable projection (VarPro) approach for solving linear and nonlinear PDEs. Neural network-based
PDE methods, especially those based on deep neural networks (DNN) [13], have flourished in the past few years; see
the review [14] and the references therein.
Variable projection (VarPro) is a classical approach for solving separable nonlinear least squares (SNLLS) problems [11,12]. These problems are separable in the sense that the unknown parameters can be separated into two sets: the linear
parameters and the nonlinear parameters. The basic idea of VarPro is to treat the linear parameters as dependent
on the nonlinear parameters, and then eliminate the linear parameters from the problem by using the linear leastsquares method. This gives rise to a reduced, but generally more complicated, nonlinear least squares problem thatinvolves only the nonlinear parameters [12]. One can then solve the reduced problem for the nonlinear parameters
by a nonlinear least squares method, typically involving a Gauss-Newton type algorithm coupled with trust region or
backtracking line search strategies [5, 2, 6, 9, 20]. Upon attaining the nonlinear parameters, one then computes thelinear parameters by the linear least squares method. Although the reduced problem is in general more complicated, thebenefits of variable projection are typically very significant. These include the reduced dimension of the parameter
space, better conditioning, and faster convergence with the reduced problem [25, 26, 12]. In some sense the idea ofInternational Conference on Scientific Computing and Machine Learning 2024 (SCML2024)

VarPro for Computational PDEs

variable projection to least squares problems can be analogized to the Schur complement in linear algebra or the staticcondensation in computational mechanics.
The VarPro algorithm was originally developed in [11], and has been improved and generalized by a number of
researchers and applied to many areas in the past few decades [15, 25, 12, 4, 22, 17, 21, 1, 3, 27, 10]. The VarPro
algorithm or its variants for training neural networks have been the subject of several studies in the literature [29, 31, 30,
26, 23, 16, 19, 18, 8]. In [26] the authors have proved that the reduced nonlinear functional of the variable projection
approach, while seemingly more complicated, leads to a better-conditioned problem and always converges faster thanthe original problem; see also [25]. The VarPro method together with the Levenberg-Marquardt algorithm is employed
for the training of two-layered neural networks in [23, 16] and compared with other related approaches. In [19, 18]
the authors extend the variable projection approach to deal with non-quadratic objective functions, and also present a
stochastic optimization method based on variable projection for training deep neural networks with attractive properties.
In this paper we focus on the variable projection approach together with ANN for solving partial differential equations.
The general strategy is as follows. We employ a feed-forward neural network to represent the field solution to the PDE,
requiring that the output layer be linear and with zero bias. We enforce the PDEs and boundary/initial conditions on aset of collocation points in the domain and on the boundary. This gives rise to a set of discrete equations about theweight/bias coefficients in the neural network. In turn, this set of equations leads to a nonlinear least squares problem
about the network coefficients.
If the boundary/initial value problem is linear, then the aforementioned nonlinear least squares problem is separable,
and the problem can be solved by VarPro. On the other hand, if the boundary/initial value problem is nonlinear, theaforementioned nonlinear least squares problem is not separable. As a result, the variable projection cannot be directly
used for nonlinear PDEs. In this case we present a combined Newton-variable projection method together with ANNsfor solving the problem. We first linearize the problem for the Newton iteration. The linearized system is linear withrespect to the updated approximation field, and it is solved by the variable projection approach together with ANNs.
Therefore, to solve nonlinear PDEs, our method involves an overall Newton iteration. Within each iteration, we use
VarPro together with ANNs to solve the linearized system to attain the updated field approximation. Upon convergenceof the Newton iteration, the neural-network coefficients contain the representation of the solution field to the original
nonlinear problem.
We present numerical examples with both linear and nonlinear PDEs to test the performance of VarPro. For smooth field
solutions, the VarPro errors decrease exponentially as the number of collocation points or the number of output-layer
coefficients increases, which is reminiscent of the spectral convergence of traditional high-order methods. We also
compare the performance of the VarPro method with that of PINN [24].

2

Variable Projection with Artificial Neural Networks for Computational PDEs

2.1

Linear PDEs

Consider the following generic linear boundary-value problem on domain Ω,
Lu = f (x), x ∈ Ω,
(1a)
Bu = g(x), x ∈ ∂Ω.
(1b)Here u(x) is the field solution to be solved for, L denotes a linear differential operator, B denotes a linear operator on
the boundary ∂Ω representing the boundary conditions, and f (x) and g(x) are prescribed non-homogeneous terms inthe domain or on the boundary.
We approximate u(x) by a feed-forward neural network, and require that the output layer be linear and contain no bias.
Let Φj (θ, x) (1 ⩽ j ⩽ M ) denote the output fields of the last hidden layer, where θ = (θ1 , . . . , θNh )T denotes thevector of weight/bias coefficients in all the hidden layers of the network. Then the NN logic leads to the following
relation,MXu(x) =
βj Φj (θ, x) = Φ(θ, x)β(2)
j=1

where Φ(θ, x) = [Φ1 (θ, x), . . . , ΦM (θ, x)] denotes the set of output fields of the last hidden layer, and β =[β1 , . . . , βM ]T is the vector of output-layer weights. Note that (θ, β) are the trainable parameters of the neuralnetwork.
We choose a set of N (N ⩾ 1) collocation points on Ω. Among them Nb (1 ⩽ Nb ⩽ N − 1) collocation points
reside on the boundary ∂Ω, and the rest of the points are from the interior of Ω. We use X to denote the set of all the
collocation points and Xb to denote the set of collocation points on ∂Ω.
2

VarPro for Computational PDEs

Table 1: Poisson equation: comparison of the maximum/rms errors and network training time by VarPro and PINN/Adam. NN
architecture [2,100,1].
collocation
points
5×5
10 × 1015 × 15
20 × 20

max-error
6.470E + 1
8.388E − 3
6.018E − 7
3.693E − 7

VarProrms-error
2.422E + 1
3.941E − 38.241E − 8
4.216E − 8train-time(secs)
1.85E + 0
5.13E + 0
2.17E + 1
2.88E + 1

PINN
rms-error
2.13E − 1
2.03E − 3
2.57E − 3
2.57E − 3

max-error
7.65E − 1
1.11E − 2
1.83E − 2
1.80E − 2train-time(secs)
3.14E + 1
1.12E + 2
1.43E + 2
1.81E + 2

Inserting the relation (2) into (1), and enforcing the equation (1a) on all the collocation points from X and the equation(1b) on all the boundary collocation points from Xb , we arrive at the following system,M
X[LΦj (θ, xp )] βj = f (xp ), 1 ⩽ p ⩽ N, where xp ∈ X,
(3a)
j=1
M
X

1 ⩽ q ⩽ Nb , where xq ∈ Xb .

[BΦj (θ, xq )] βj = g(xq ),

(3b)

j=1

This is a system of (N + Nb ) algebraic equations about the trainable parameters (θ, β), with (Nh + M ) unknowns.
We seek a least squares solution for (θ, β) to the system (3). This system is linear with respect to β, and nonlinear withrespect to θ. This leads to a separable nonlinear least squares problem. To make the formulation more compact, we
re-write the system (3) into a matrix form,
H(θ)β = S,
(4)
(N +Nb )×M(N +Nb )where H(θ) ∈ M
and S ∈ R
. For any given θ, the least squares solution for the linear parameters β+
to this system is given by β = [H(θ)] S, where the superscript in H+ denotes the Moore-Penrose pseudo-inverse of
H. Define the residual of the system (4) byr(θ) = H(θ)β − S = H(θ)H+ (θ)S − S.
(5)
We compute the optimal nonlinear parameters θopt by minimizing the Euclidean norm of r,
1
1
θopt = arg min ∥r(θ)∥2 = arg min ∥H(θ)H+ (θ)S − S∥2 .
(6)
2
2
θ
θ
After θopt is obtained, we can compute the optimal linear parameters βopt by solving equation (4) using the linear leastsquares method. The problem represented by (6) is a nonlinear least squares problem about θ only, where the linear
parameter β has been eliminated. We solve this problem by a Gauss-Newton algorithm combined with a trust region
strategy [6, 9, 7, 28].
2.2

Nonlinear PDEs

Consider the following nonlinear boundary value problem,
Lu + F (u) = f (x),(7a)
Bu + G(u) = g(x), on ∂Ω,(7b)where F (u) and G(u) are nonlinear operators on the solution field u(x) and also possibly on its derivatives, and L, B,
f and g have the same meanings as in the equations (1a)–(1b).
To exploit VarPro, we first linearize the system (7a)–(7b) with the Newton’s method. Let uk denote the approximationof the solution at the k-th Newton iteration. We linearize this system as follows,
Luk+1 + F (uk ) + F ′ (uk ) uk+1 − uk = f (x),
(8a)

k+1
k
′ k
k+1
k
Bu
+ G(u ) + G (u ) u− u = g(x), on ∂Ω,(8b)′′where F (u) and G (u) denote the derivatives with respect to u. We further re-write the linearized system into,
Luk+1 + F ′ (uk )uk+1 = f (x) − F (uk ) + F ′ (uk )uk ,
Bu

k+1

′

k

k+1

+ G (u )u

k

′

k

k

= g(x) − G(u ) + G (u )u ,k

(9a)on ∂Ω.

(9b)
k+1

Given u , this system represents a linear boundary value problem about the updated approximation field u. Therefore,
the VarPro/ANN algorithm from Section 2.1 can be used to solve this linearized system (9a)–(9b) for uk+1 . Upon
convergence of the Newton iteration, the solution to the original nonlinear system (7a)–(7b) will be obtained andrepresented by the neural-network coefficients.
3

VarPro for Computational PDEs

3

Numerical Examples
10

3.1

Poisson Equation

1

10

-1

10

-3

max error
rms error

We consider the 2D Poisson equation ∇2 u =f (x, y) with Dirichlet
BCs on the unit squaredomain Ω = [0, 1] ×
Collocation points per direction
(a)
(b)
(c)
[0, 1], with an analytic
solution as shown in Figure 1: Poisson equation: (a) exact solution. (b) Pointwise error of the VarPro solution. (c)Fig. 1(a). Fig. 1(b) de- Maximum/rms errors of VarPro versus the number of collocation points per direction. NN architecture
picts the point-wise error [2, 20, 100, 1], cos activation function.
distribution of the VarPro solution obtained on an NN architecture [2, 20, 100, 1] with the cosine activation function,
showing a maximum error on the order 10−9 in the domain. Fig. 1(c) shows the the VarPro maximum/rms errors versus
the collocation points, indicating an exponential convergence. Table 1 compares the maximum/rms errors and the
network training time between VarPro and PINN, showing that VarPro is much more accurate and computationallymore efficient.
10

-5

10-7

10

10

3.2

-9

-11

5

10

15

20

25

30

Viscous Burgers Equation

We next consider theviscous Burgers equationon the domain (x, t) ∈[−1, 1] × [0, 1.05],
∂u
∂u
1 ∂2u
∂t + u ∂x = 100π ∂x2 ,with homogeneous BC
and the initial conditionxu(x, 0) = − sin(πx).
(a)
(b)
(c)
This problem has asolution in which a sharp Figure 2: Burgers equation: distributions of (a) VarPro solution, and (b) its point-wise error. (c) Comgradient (akin to a jump parison of profiles between the VarPro solution and exact solution at t = 1. Domain decomposition
discontinuity) develops and block time marching are used.
over time. Figs. 2(a,b) show distributions of the VarPro solution and its point-wise absolute error in the domain, with amaximum error on the order 10−5 . Fig. 2(c) compares the solution profiles between the VarPro solution and the exact
solution at t = 1, illustrating the sharp gradient in the domain. Block time marching and domain decomposition (with 6
sub-domains) have been used with VarPro for this problem, with a local NN architecture [2, 250, 1] on each sub-domain.
VarPro result is highly accurate.
1

VarPro solution
Exact solution

Profiles

0.5

0

-0.5

-1
-1

4

-0.5

0

0.5

1

Concluding Remarks

We have presented a variable projection-based method together with artificial neural networks for solving linear andnonlinear partial differential equations. The basic idea of variable projection (VarPro) is to distinguish the linear
parameters from the nonlinear parameters, and then eliminate the linear parameters to attain a reduced formulation of
the problem. One can then solve the reduced problem for the nonlinear parameters first, and then compute the linear
parameters by the linear least squares method afterwards. For smooth field solutions, the errors of the VarPro method
decrease exponentially or nearly exponentially with increasing number of collocation points or with increasing numberof output-layer coefficients. The test results show that the VarPro method is highly accurate. Even with a fairly small
number of nodes in the neural network, or with a fairly small set of collocation points, the VarPro method can produce
very accurate simulation results.

Acknowledgments
This work was partially supported by US National Science Foundation (DMS-2012415).
4

VarPro for Computational PDEs

References[1] T. Askham and J.N. Kutz. Variable projection methods for an optimized dynamic mode decomposition. SIAM J.
Applied Dynamical Systems, 17:380–416, 2018.
[2] A. Bjorck. Numerical Methods in Matrix Computations. Springer, 2015.
[3] G.-Y. Chen, M. Gan, C.L.P. Chen, and H.-X. Li. A regularized variable projection algorithm for separable
nonlinear least-squares problems. IEEE Transactions on Automatic Control, 64:526–537, 2019.
[4] J. Chung, E. Haber, and J. Nagy. Numerical methods for coupled super-resolution. Inverse Problems, 22:1261–
1272, 2006.
[5] J.E. Dennis and R.B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations.
SIAM, 1996.
[6] S. Dong and Z. Li. Local extreme learning machines and domain decomposition for solving linear and nonlinear
partial differential equations. Computer Methods in Applied Mechanics and Engineering, 387:114129, 2021. (alsoarXiv:2012.02895).
[7] S. Dong and Y. Wang. A method for computing inverse parametric pde problems with random-weight neuralnetworks. Journal of Computational Physics, 489:112263, 2023. (also arXiv:2210.04338).
[8] S. Dong and J. Yang. Numerical approximation of partial differential equations by a variable projection methodwith artificial neural networks. Computer Methods in Applied Mechanics and Engineering, 398:115284, 2022.
(also arXiv:2201.09989).
[9] S. Dong and J. Yang. On computing the hyperparameter of extreme learning machines: algorithms and applications
to computational PDEs, and comparison with classical and high-order finite elements. Journal of Computational
Physics, 463:111290, 2022. (also arXiv:2110.14121).
[10] N.B. Erichson, P. Zheng, K. Manohar, J.N. Kutz S.L. Brunton, and A.Y. Aravkin. Sparse principal component
analysis via variable projection. SIAM J. Appl. Math., 80:977–1002, 2020.
[11] G.H. Golub and V. Pereyra. The differentiation of pseudo-inverse and nonlinear least squares problems whose
variables separate. SIAM J. Numer. Anal., 10:413–432, 1973.
[12] G.H. Golub and V. Pereyra. Separable nonlinear least squares: the variable projection method and its applications.
Inverse Problems, 19:R1–R26, 2003.
[13] I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. The MIT Press, 2016.
[14] G.E. Karniadakis, G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. Physics-informed machine learning.
Nature Reviews Physics, 3:422–440, 2021.
[15] L. Kaufman. A variable projection method for solving separable nonlinear least squares problems. BIT, 15:49–57,
1975.
[16] C.-T. Kim and J.-J. Lee. Training two-layered feedforward networks with variable projection method. IEEE
Transactions on Neural Networks, 19:371–375, 2008.
[17] K.M. Mullen and I.H.M. van Stokkum. The variable projection algorithm in time-resolved spectroscopy, microscopy and mass spectrometry applications. Numer. Algor., 51:319–340, 2009.
[18] E. Newman, J. Chung, M. Chung, and L. Ruthotto. SlimTrain – a stochastic approximation method for training
separable deep neural networks. arXiv:2109.14002, 2021.
[19] E. Newman, L. Ruthotto, J. Hart, and B. van Bloemen Waanders. Train like a (Var)Pro: Efficient training ofneural networks with variable projection. arXiv:2007.13171, 2020.
[20] N. Ni and S. Dong. Numerical computation of partial differential equations by hidden-layer concatenated extreme
learning machine. Journal of Scientific Computing, 95:35, 2023. (also arXiv:2204.11375).
[21] D.P. O’Leary and B.W. Rust. Variable projection for nonlinear least squares problems. Comput. Optim. Appl.,
54:579–593, 2013.
[22] M.R. Osborne. Separable least squares, variable projection, and the gauss-newton algorithm. Electronic Transactions on Numerical Analysis, 28:1–15, 2007.
[23] V. Pereyra, G. Scherer, and F. Wong. Variable projections neural network training. Mathematics and Computersin Simulation, 73:231–243, 2006.
[24] M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: a deep learning framework forsolving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational
Physics, 378:686–707, 2019.
5

VarPro for Computational PDEs[25] A. Ruhe and P.A. Wedin. Algorithms for separable nonlinear least squares problems. SIAM Review, 22:318–337,
1980.
[26] J. Sjoberg and M. Viberg. Separable nonlinear least squares minimization - possible improvements for neural net
fitting. Neural Networks for Signal Processing VII. Proceedings of IEEE Signal Processing Workshop, 1997.
[27] X. Song, W. Xu, K. Hayami, and N. Zheng. Secant variable projection method for solving nonnegative separable
least squares problems. Numerical Algorithms, 85:737–761, 2020.
[28] Y. Wang and S. Dong. An extreme learning machine-based method for computational pdes in higher dimensions.
Computer Methods in Applied Mechanics and Engineering, 418:116578, 2024.
[29] K. Weigl and M. Berthod. Neural networks as dynamical bases in function space. Report No 2124, INRIA,
Sophis-Antipolis, France, 1993. URL: https://hal.inria.fr/inria-00074548/document.
[30] K. Weigl and M. Berthod. Projection learning: alternative approach to the computation of the projection. Proc.
European Symp. on Artificial Neural Networks, Brussels, Belgium, pages 19–24, 1994.
[31] K. Weigl, G. Giraudon, and M. Berthod. Application of projection learning to the detection of urban areas in SPOT
satellite images. Report No 2143, INRIA, Sophia-Antipolis, France, 1993. URL: https://hal.inria.fr/inria-00074529.

6