Deep Learning for Subspace Regression
Authors: Vladimir Fanaskov
Vladislav Trifonov
Alexander Rudikov
Ekaterina Muravleva
Ivan Oseledets
Venue: ICLR 2026 Poster (2026)
Keywords: grassmannian
regression
subspace regression
supervised learning
ROM
POD
optimal control
balanced truncation
parametric PDEs
eigenproblems
deflated conjugate gadient
coarse grid correction
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: Regression on the Grassmann manifold works well with deep learning.

Abstract: It is often possible to perform reduced order modelling by specifying linear subspace which accurately captures the dynamics of the system. This approach becomes especially appealing when linear subspace explicitly depends on parameters of the problem. A practical way to apply such a scheme is to compute subspaces for a selected set of parameters in the computationally demanding offline stage and in the online stage approximate subspace for unknown parameters by interpolation. For realistic problems the space of parameters is high dimensional, which renders classical interpolation strategies infeasible or unreliable. We propose to relax the interpolation problem to regression, introduce several loss functions suitable for subspace data, and use a neural network as an approximation to high-dimensional target function. To further simplify a learning problem we introduce redundancy: in place of predicting subspace of a given dimension we predict larger subspace. We show theoretically that this strategy decreases the complexity of the mapping for elliptic eigenproblems with constant coefficients and makes the mapping smoother for general smooth function on the Grassmann manifold. Empirical results also show that accuracy significantly improves when larger-than-required subspaces are predicted. With the set of numerical illustrations we demonstrate that subspace regression can be useful for a range of tasks including parametric eigenproblems, deflation techniques, relaxation methods, optimal control and solution of parametric partial differential equations.

Published as a conference paper at ICLR 2026
DEEP LEARNING FOR SUBSPACE REGRESSION
Vladimir Fanaskov
Vladislav Trifonov A14Science, AIC
Alexander Rudikoy AXXX, INM
AXXX
fanaskov.
wwadimirr@mailil.......i.iiaaa.a fanaskov..laaddiiiiooo com
Ekaterina Muravleva AI4Science, AIC
Ivan Oseledets AXXX, INM
ABSTRACT
It is often possible to perform reduced order modelling by specifying linear sub- space which accurately captures the dynamics of the system. This approach be- comes especially appealing when linear subspace explicitly depends on parame- ters of the problem. practical way to apply such a scheme to compute sub- spaces for selected set of parameters in the computationally demanding offline stage and in the online stage approximate subspace for unknown parameters by interpolation. For realistic problems the space of parameters is high dimensional, which renders classical interpolation strategies infeasible or unreliable. We pro- pose to relax the interpolation problem regression, introduce several loss fune- tions suitable for subspace data, and use a neural network as an approximation nigh-dimensional target function. To further simplify learning problem we introduce redundancy: in place of predicting subspace of a given dimension we predict larger subspace. We show theoretically that this strategy decreases the complexity of the mapping for elliptic cigenproblems with constant coefficients and makes the mapping smoother for general smooth function on the Grassmann manifold. Empirical results also show that accuracy significantly improves when larger-than-requrre subspaces are predicted. With the set numerical illustra- tions we demonstrate that subspace regression can be useful for range of tasks in- cluding parametric eigenproblems, deflation techniques, relaxation methods, op- timal control and solution of parametric partial differential equations.
INTRODUCTION
The goal of reduced order modelling (ROM) to identify uninformative degrees of freedom and discard them (Bai al. 2005). The result is simplified system that is casier to analyse and simu- late. This program is computationally demanding and only justified in the setting when many related problems repeatedly solved and it is possible use information from encountered problems to build reduced model for the problems to come. Typical examples are parametric models for par- tial and ordinary differential equations (PDEs and ODEs), and usual applications are optimisation, sensitivity analysis, uncertainty quantification and control.
As an illustration consider proper orthogonal decomposition (POD) for time-dependent PDEs Volk- wein 2013). (Hesthaven et al. 2022). To apply POD, one computes solutions for representative set of parameters and builds reduced basis for spatial variables by the best low-rank approximation. When new parameters arrive, computed basis is used to discretise PDE that is solved at reduced cost. In global POD this basis is the same for all incoming parameters and in local POD basis ex- parameters. As one may expect, local POD is more expressive than global
plicitly depends POD, but can be more challenging to arrange.
POD is an example of the general class of techniques where linear subspace parametrises useful degrees of freedom. While nonlinear ROM techniques exist, linear methods are better understood theoretically, easier apply in practice and provide sufficiently well approximation, especially when local versions are available (Franco et al., 2024). In this setting the main challenge to construct reliable approximation the function that maps new parameters to linear subspaces.

---

Published as a conference paper at ICLR 2026
We analyse this problem under the following assumptions: (i) the set of parameters of interest is specified in form probability distribution, (ii) it known how to compute good or optimal linear subspace for each parameter, (iii) the numerically stable method to construct reduced problem from basis is available. In short, we consider regression on grassmannian. We approach the regression problem by specifying loss function and using neural networks parametric models accommo- date high-dimensional parameter spaces pervasive in practical problems.
More specifically, our main contributions are:
1. Mathematical formulation of subspace regression problem and examples of applications, including eigenspace approximation, local POD, learning basis for deflation and two-grid method, approximating salanced-truncation basis for optimal control problems.
2. Several loss functions, suitable for neural network training, including the stochastic one that scales well with the increase of subspace size.
3. Embedding technique: a strategy to leam larger subspace containing the target one. Em- pirically, this technique significantly improves accuracy for subspace learning.
4. Two theoretical justification embedding technique: derivative of smooth function on Grassmann manifold can be reduced by embedding; complexity of the subspace regression problem for elliptic cigenproblem with constant coefficients.
5. Empirical evaluation of proposed techniques on diverse problems including com- parisons with neural surrogates, kernel methods and classical interpolation in normal coor- dinates.
SUBSPACE REGRESSION
In this section we formulate precisely what we mean by subspace regression and describe several applications.
2.1 DEFINITION OF SUBSPACE REGRESSION PROBLEM
In lincar space R" we define k-dimensional subspace S(W) = {Wa, a E R²} full rank matrix W' E R R Matrices W1 and W2 represent the same subspace if there is an invert- ible matrix G such that W1 WzG. The equivalence class of such matrices is denoted by [W]. The set of all k-diiensional subspace of n--dimensional space Gr(k, is known as Grassmann manifold or grassmanian (Ciaramella et al.) 2025). (Bendokat et al.) 2024).
by specifying tall
Let V IR + Gr(k, 12) be a function that maps the space of parameters r E R² to the subset of grassmanian represented as the set of tall full rank matrices V(r). We assume that parameters are sampled from distribution r Pr and that dataset 'D = {(r1, 1/1,,.. ...... (Fms Vm)} of m i.i.d. samples is available. For a given parametric model Yo : IP - Gr(r, n),r 2 we want to identify parameters 0* such that Yo- (r) approximates (rr We formulate this task as optimisation problem
8* = arg min /
(LY(r)) ~ arg min prupr
E
Σ 1=1
(1)
Loss function for subspace regression problem is assumed to have two properties:
L(A, L(Ā, B) for arbitrary E [A], B E [B]: L(A, B) and L(A, B) iff S(B) C S(A).
(2)
In Section 3 we provide explicit expression for loss functions with these properties.
Aside from unusual invariance requirement (2). optimisation problem (T) standard machine learning formulation of regression problems which can be solved with stochastic optimisation for arbitrary model (r(r) that admits efficient evaluation of gradients.
Note, that we allow Yo to have more columns than target V. In this context approximation is understood in terms of subspace inclusion S (V) C S (Wo). As we explain later, redundancy introduced this way significantly improve accuracy.
2

---

Published as a conference paper at ICLR 2026
2.2 EXAMPLES OF SUBSPACE REGRESSION PROBLEM
Approximate eigenspaces. Consider eigenproblem for Schrödinger equation
-Av(z) + (x)) Ex)x), || < 00:
(3)
where Uxx) is potential energy and E' is energy of the system. One way to find eigenpairs is to ap- proximate (zx) by a finite series (xx)= = Ek (2) and enforce Petrov-Galerkin condition that residual is orthogonal to all (x(x)). Continuous problem (3) reduces to cigenproblem for Hermitian matrix and can be solved in O (K*) operations (Trefethen & Bau 2022). For eigenproblems we are typically interested only in extremal eigenspaces corresponding to either smallest largest eigen- values. In this case it is desirable to select small number of basis functions Φ2 (zx) that approximate sufficiently well the subspace of interest. When eigenproblem (3) is solved repeatedly for many potential functions (z) this lead us to subspace regression problem ! used to approximate the
mapping (xx) + [(x) [(x)= = Σ (nd)(r), a; E C That is, we wish to predict subspace spanned by first K eigevectors. When this mapping learned from a set of examples, eigenproblems for unobserved potentials U can be solved efficiently, since www-dimesiiona candidate subspace for eigenfunctions is available.
Intrusive POD for time-dependent PDEs. As an example of time-dependent PDE we use Burgers equation
tu(a,t) t) 1(x,1)) Øt
tu(x, 1) 0
Ox Ox (v(x) (((x) (ute))))))) u(0, =000,,)==(11))) 0, u(x; 0) un(). (4)
)u(z,()
Ox
One starts with spatial discretisation which reduces equation to the set of ODEs and de- fine inner product (::)w for discretised u(t) that approximates Le inner product. For this set of ODEs the reduced degrees of freedom v are defined as solution to optimisation problems min S dt [lu(t) - (w,, u(t)) Willit subject to (ih; 45)} = 0, < i, (wh,, v/a)
When dis- cretised, this scheme lead to optimal basis computed with SVD from snapshot matrix (Volkwein) 2013). This basis can only be computed when equation (4) is integrated, so POD is justified only in situation when many related problems are solved. We apply subspace regression with POD to learn the function that maps PDE data to the subspace formed by reduced basis {1),... th for some small k. Notably, this allows to apply local POD to high high-dimenssionl parametric problems.
Coarse grid correction for iterative methods. Consider stationary diffusion equation with Dirich- let boundary conditions
(5) When equation 55 discretised with finite difference or finite element method, reduces to linear problem Aø = T. where A is large sparse matrix and o, f are discretised solution and right-hand side of (5) To exploit sparsity of A one can solve linear equation with relaxation method. General relaxation method split matrix additively A = D + C, where D is regular with known inverse [Saad, 2003). Given the split, if iteration scheme n==+1 = x. + D-1 (b Ax") is convergent,
dii k(x) grad d(x) = f(x), x ET, ((z)aar = 0.
steady state is exact solution x A b Convergence is linear and its rate is defined by spectral radius of error propagation matrix I - D-1A, To improve convergence rate, one can augment iterative method with coarse-grid correction (Trottenberg al. 2001). This techniques allows one to remove influence of leading subspace formed by colums of matrix V of T D- A by solving small reduced linear system for error equation V AVe r, where and are error and residual in the subspace S(V). Naturally, subspace regression (I) for this problem approximates the mapping A + S(V or k(x)) ~ S(V) for linear systems resulting from equation
Deflation for conjugate gradient. Krylov subspace methods provide a more systematic way to solve large sparse linear systems (Saad) 2003). For linear system with symmetric positive definite matrix A resulting from discretisation of equation (5). the method of choice is conjugate gradient (CG) (Hestenes et al. 1952). Similar to other Krylov methods, on step r, CG identify optimal solution within Krylov subspace K, span (b, Ab, Ab. Ab-b),, where span refers to the subspace formed by linear combinations of vectors in the set. Since powers of A are involved, the most readily available vectors are from the subspaces with large eigenvalues (Saad 2011). To improve
Suitable discretisation of parametrisation of both & and U is assumed.
3

---

Published as a conference paper at ICLR 2026
convergence of method it is reasonably include eigenspaces V with small eigenvalues to the approximation space K.. The resulting method is deflated CG and the approximation space is K,, U S(V) (Saad et al. 2000). In analogy with previous example, the subspace regression problem considered approximates A -> S(V) or k(x) - S(V), but this time spans eigenspaces with small eigenvalues of matrix A.
dalanced-truncatiin for linear-quadratic control. Suppose we want to solve the following linear quadratic control problem
j((t) Ay(t) + Bu(t), 2(t) = Cy(t),
J J= f dt ((=(t)) Q2(t) + (u(t)) Ru(t)) + (=(T)) Mz(T),
where vtt is state variable, u(t) is control, z(t) is observable, A, A, Q, R, M matrices of appropriate sizes, R is symmetric positive definite, Q and M are symmetric positive semidefinite. The goal is to find control signal u((), I [0, 7] that minimises cost function J.
In the situation when number of state variables y)1) is large, one may want to apply ROM to compute optimal control at a reduced cost. An established way to do that is balanced truncation (Moore 2003). Roughly speaking, balanced truncation compute special coordinate system y(t) = Ty()) that discounts variables that are both unobservable and uncontrollable, so only a few first columns of matrix T can be used to accurately model (6). This is done by finding coordinate system that simultaneously diagonalises observability gramian G, and controlability gramian G, defined as solutions of Lyapunov equations A W. + W.A W C C 0, AW. + W. + BB 0 (Moore) 2003), (Volkwein) 2013). In this case the goal of subspace regression ① is to approximate the mapping A, B, (T). where T is tall matrix assembled from first few columns of T.
3
THEORETICAL RESULTS
We proceed by characterising loss functions, introducing the subspace embedding technique and providing its theoretical justification.
3.1 LOSS FUNCTIONS
Requirements 2 that allow loss function to work with Gr(k, n) data enforce right GL.(k) invariance. As consequence all loss functions introduced below are all based on orthogonal projectors.
Theorem 1. Let A E Rank R=xp,p S tall full rank matrices.
1. Loss function L) = p- 1- satisfies requirements 2 where A
QARA. B QgRa are reduced QR decompositions of A and B. IF TF is Frobenius norm
2. Let z & be random variable with zero mean and identity covariance matrix. Loss functions L2(A, B;2) min, Adu-Qp=llp does not depend on the choice of A from
[A], where B = QBRp 3 is QR decomposition.
3. On average L2 equals L,, i.e., E. [L, B;2)] ().
Proof. Appendix
Loss L1 is essentially the same as the difference of orthogonal projectors. Note, that (A, with equality reached if and only if matrices A and B share the same columns space, since in this case [[Q distance is replaced with error of least squares problem; (ii) to remove second projector, stochastic Hutchinson trace estimation used. Reformulation with least square problem allows one use normal equation, and various tools from randomised numerical linear algebra, e.g., randomised
p. Loss Ta introduces two modifications: (i) projector in Riemannian
Reduced QR decomposition of tall full rank matrix A Baxk is factorisation A RA. where QA x has orthonormal columns, R. TR is upper triangular with nonzero elements the diagonal.
4

---

Published as a conference paper at ICLR 2026
preconditioned Cholesky- -RR (Garrison & Ipsen 2024). blendenpik solver (Avron et al. 2010), and sketching [Woodruff et al. 2014). We will see in Section 4 that loss function L2(A, B) based on normal equation scales better than L1 (A, B) with the increase of subspace size.
3.2 SUBSPACE EMBEDDING
In the definition of subspace regression problem we allow to approximate function R' - Gr(k, n) by function E' - Gr(r,n) where It is justified because of two unique properties of regression and interpolation inclusion of vector subspaces well-defined; (ii) subspace, predicted by regression model or inter- polated by standard techniques. is used to construct a reduced model. From the latter property one may expect similar or improved accuracy when the predicted subspace from Gr(r, n) contains target subspace from Gr( n.).
We call this strategy subspace embedding. grassmanian: (i)
We will show empirically in Section 4hhtt that subspace embedding significantly improves accuracy and generalisation gap. Here argue that prediction of larger-than teerr-tthh--nedddd subspaces align well with inductive bias of neural networks known as f-principle or spectral bias (Xu et al. 2019). F-principle an observation that neural networks tend to learn smoothed versions of the target functions. As we show below, embedding techniques may improve smoothness of learned function.
Theorem 2. Let V R - Gr(k, n) be continuously differentiable on t E [0, V(t) V(t) = 4s- It is always possible construct piecewise continuous function W I + Gr(r, r > k. ( (t) = I, such that 2/2 ([) - (t))
- r-k is arbitrary small and
w $ V(t where inequality strict for all points where | #0.
Proof. Appendix See Appendix C for subspace embedding example.
Theorem implies that, by increasing the subspace size, one can always approximate continuously differentiable functions arbitrarily well and simultaneously reduce its derivative. F-principle sug- gests that the latter property makes learning easier for neural networks.
3.3 COMPLEXITY OF PARAMETRIC EIGENPROBLEM
To illustrate difficulties may encounter and further justify embedding technique consider complexity of subspace regression problem for parametric elliptic eigenproblem with constant co- efficient
D
82
Ph...kk in (1) (....... ID)
Σa
a
Xin....in i (x (#1,.... ................................ TD),
Oz?
where a; > 0, 14 E [0, 1] and Dirichlet boundary conditions are assumed.
7 general cigenfunction is (F = I sin (TiJH,)) and the set of eigenfunctions does not depend on coefficients 04- Observe that Xip....ip = Σ aj (miz)² so coefficients a; define the order of eigenvectors. Below we formally characterise mapping from coefficients to k--th eigenvector and eigenspace.
For problem
Theorem 3. Suppose eigenvectors of 7 are ordered according to the increase of eigenvalues. Let du be an eigenvector on position k, let Vv be an eigenspace spanned by vectors on positions up to k. Consider mappings Fr. aa, + dn and G1. "l,...,ap + Vv-
1. Fg, GE are piecewise constant functions that map real numbers to elements of sets SF. SG.- Sets SF..SG, are finite with # (ks, D), #G (k, D) distinct elements.
2. Let Wi, > be a subspace obtained by union of Vi for distinct .p. Number of vectors in W, is # (1, D) +
1 3. # (k. D) (D-1)1 (log k)) D-1 for fixed D and large k.
4. # r, (k: D) k.Dloga & for fixed k and large D.
5

---

Published as a conference paper at ICLR 2026
5. #o. () 2 (D-1) 1D-1IkD-1 1 for fixed D and large k.
kD
6. #o (k, D) 2 -ύ- D* 1 for fixed is and large D.
Where ~ is asymptotic expansion and 2,5 are lower and upper bound on leading asymptotic.
Proof. Appendix
D See Appendix for examples.
Theorem suggests that for problem mappings from coefficient to k-th eigenvector de or sub- space Vi are piecewise constant functions with rapidly growing number of constant regions when either & or D increases. The complexity (the number of regions) of the subspace prediction prob- lem exceeds the complexity of k-th eigenvector prediction. However, results also suggest that the number of unique eigenvectors within V& grows at the same rate as the number of eigenvectors on position k. This means, the large number of distinct regions in G' comes from large number of possible combinations of an asymptotically small number of vectors. Given that, the complexity of mapping from coefficients to subspace decreases, if one predicts subspace of larger dimension Vi 2 Vg- For example, one is willing to predict subspace of dimension # p (k, D) in place of Vx of dimension kt the mapping CI....,ap + V2 2 V₂ may be chosen to have constant value.
4 NUMERICAL EXPERIMENTS
We present several numerical experiments to corroborate our theoretical findings. The discussion of control problems appears in Appendix [] All numerical results are reported as single-run metrics without explicit егrоr bars. To study sensitivity to train-test split we perform several dedicated ex- periments. Results suggest that variance is low and does not affect main conclusions. The details available in Appendix [] For most problems predicted quantity e and ground truth value e it reads 100% x ||v ollz ||0||2- For eigenvalue prob- lems the numerator is Z2 adjusted explained below. For iterative methods report convergence plots for relative erro, without additional factor 100%.
report relative error measured in percents. For
4.1 EIGENSPACE PREDICTION
We considered several eigenvalue problems: (i) D eigemproblem with Schrödinger operator (3)
with parametric family of expanded Morse oscillator (Le Roy et al. 2006). (ii) D = 2 Schrödinger
operator (3) with parametric family of two expanded Morse oscillators (Carpenter et al. 2018).
(ii) D = 2 datasets, kn li2 and kn k2, for elliptic eigenproblem (5) (left-hand side of
the equation) with contrast coefficient sampled from gaussian random field, (iv) D = 3 dataset for elliptic eigenproblem with diffusion coefficient k kg- In all experiments we used FFNO architecture (Tran et al., (202T) a modification of FNO (Li et al. (020), and performed extensive hyperparameter grid search. Details on dataset generation and training protocol are available in Appendix F To contextualise subspace regression we provide results for two baselines.
Regression with Zg adjusted le loss. Eigenvectors are defined up sign, so in place of subspace losses specified in Theorem [ one can try to directly predict eigenvectors with Z, adjusted l2 loss 1z, (v, = min ||u wllp-
Interpolation in Riemannian normal coordinate system. A standard technique of manifold in- terpolation applied to grassmannian (Amsallem) 2010), (Ciaramella et al.) (025). (Zimmermann] 2019). For given query, k closest points are selected from the training set. One point supplies common tangent space, i.c., it is used to compute logarithms for the remaining points. Since log- arithms lay in the same tangent space they can be interpolated with any techniques desirable (we use RKHS (Bishop Nasrabadil 2006)). After interpolation of logarithms, the exponential map is computed.
Additional results are available in Appendix F Here we highlight several interesting trends.
Subspace losses unsuccessful without subspace embedding technique. Figure Ib contains results of learning subspace spanned by first 10 eigenvectors for D = 2 on grid 100 x T00 elliptic eigenproblem with Tzz( A, B;z) (loss L, (A, B) leads to the same accuracy). Neural network predicts

---

Published as a conference paper at ICLR 2026
18
L,A,, 8)
test train
90
training time. s/epoch
15
L2(A,
25
75
12
relative error,
reltive erron
35
15
3
10
20
30
40
50
60
20 30
40
6
subspace dimension
subspace dimension
elgenvector
(a) training time
(b) subspace embedding
(c) Ze adjusted l2
Figure Selected results for eigenspace prediction: (a) Comparison of training time for losses La (A, B), Tz2 (A B) from Theorem On D 2 grid N₂ = N₂ 100 we observe (A, =) scales better with dimension size; (b) Illsttationn subspace embedding technique from Section 3.2 for D elliptic eigenproblem, prediction of first 10 eigenvectors. Prediction of larger subspace manifestly improves accuracy and reduces generalisation gap; (c) Relative error for individual eigen- vector predictions for the same problem in (b) but trained with Zy Edjausseted la loss. Similarly to results of Theorem we observe steep increase problem complexity with eigenvector number. See Section 4.1 for details.
subspace of sizes 10, 20, 30, 40 according to subspace embedding strategy Section 3.2 Results demonstrate subspace embedding is efficient in decreasing test error from 30% for subspace of dimension 10, to test error 2% for subspace of dimension 40 (0.4% of the total number of degrees of freedom). is less clear from Figure lb] but the generalisation gap also systematically improves, suggesting that complexity of the problem decreases. Similar conclusions are valid for D = QM and eigenproblems.
Classical interpolation is not competitive. In Ta- ble [ we gather results (best for each method) for one D and two D QM problems. Clas- sical interpolation is reasonably accurate only on but Lz, loss still
Table 1: Relative errors for QM problems.
dataset
interp.
Lgz
L1 L1 B)
the simplest problem in D results in better accuracy. The reason is likely clas- sical interpolation struggles in high-dimensional subspaces because observations are too sparse to naively approximate tangent space in the region of interest by finding nearest neighbours. On QM datasets accuracy of subspace regression consistently better than for other approaches.
D
4.69%
2.33%
0.09%
D = 2,
31.9%
19.52%
0.65%
D 2,6
92.64%
48.56%
15.58%
Loss without QR scales much better for larger subspace sizes. In Figure la we demonstrate wall clock training time for Li (A, B) and L2(A, B; (least squares problem is solved with normal equation) per epoch on the same hardware for the same FFNO architectures, For small subspace sizes the training time the methods are roughly on par, but with the increase of subspace size QR starts to drastically slow down training with L, loss.
Training with loss L2, is reasonable only for several first eigenvectors. In Figure Ic we present results for learning individual eigenvectors (a separate network is trained for each eigenvector) for D elliptic eigenproblem. Train error reasonably small for all eigenvectors, which imply neu- ral networks successfully approximate them. Rapid growth of the test error with eigenvector number indicates the increase of problem complexity in agreement with Theorem 3
Neural networks trained with subspace embedding technique learn smoother maps. Theorem suggests it is possible to decrease derivative by embedding of geodesics into a larger space. This provides only circumstantial evidence that the may happen when neural networks trained with subspace embedding technique. In Appendix F we gather empirical results that support such conclusion. The results based several "smoothness indicators": the error of linear model,

---

Published as a conference paper at ICLR 2026
Frobenius norm of derivative, and mean cosine of angles between subspaces at nearby points. We refer interested readers to Appendix F.5
Subspace regression speed-up classical iterative eigensolvers. As example of hybrid ap- proach we consider combination of subspace regression and LOBPCG (Knyazev) 200T). LOBPCG classical iterative matrix-free eigensolver that can approximate extremal eigenspaces. To apply in combination wit hsubspace regression we use trained neural network to predict subspace, and initialise LOBPCG with solution of reduced eigenproblem. Note, that the cost of such initialisa- tion is negligible small comparing to the full cost of LOBPCG iterations. We observe 2 to 3 times faster converges and orders lower relative eтrоr on average when subspace regression is used for initialisation. More details are available in Appendix F.6
Loss L2(A, 2) may become unstable. From results summarised in Figure la one can assume that L A, B; 2)
Table 2: Accuracy for D = elliptic eigenproblem.
always preferable. Results for D = 3 problem (el- liptic eigenproblem, grid 30 x 30 x 300 pre- diction of first 3 eigenvectors) summarised
Nanh
LI(A,
L (A, B;
6
24.77%
31.46%
28.28%
in Table 2 reveal a more nuanced picture. Loss L2 L2,A, =) clearly performs worse than L1 (A, B) and even fails for subspace size 24. The reason for that is numerical instability of solvers based on the normal equation. To stabilize L2(A, 2) apply Cholesky-QR2 (Yamamoto et al. 2015). The results for stabilised loss shows that accuracy becomes comparable to L1 B) and even slightly better for larger subspace dimensions.
13.69%
17.12%
15.88%
24
9.71%
9.49%
48
7.54%
16.3%
7.4%
4.2 PARAMETRIC PDE PROBLEMS
We considered two PDEs: (i) == 1 + 1 viscous Burgers equation, related to benchmark from
al. 2020); (ii) = 2 elliptic problems ( Our main operator is FFNO and the solutions strategy
we use is classical intrusive P For datasets description and training details see Appendix G We compare subspace regression with several methods.
Regression with FFNO. We apply FFNO, extension of Fourier Neural Operator, to parametric PDEs in standard way similar to Tran et al. 2021
Regression with DeepONet. Classical architecture based on the universal approximation theorem of operators (Lu et al. 2019). DeepONet can be understood as end-to-end training of non-intrusive POD with basis functions parametrised by implicit neural representation or physics-informee neural networks (Sitzmann et al. 2020), (Lagaris et al., 1998), (Raissi et al. 2019).
Intrusive POD with DeepONet/FFNO basis. When DeepONet is trained, learned spatial spa- tiotemporal basis functions can be used to extract basis (by method directly related to POD) suitable for spectral methods or intrusive POD (Meuris et al., 202I), (Meuris al. (202)). As suggested in the discussion section of (Meuris et al., 2021). the same can be done with FNO.
Deep POD. Projector- based loss is used directly extract basis from available snapshot matrices or steady-state solutions (Franco et al.). In the referenced publication this approach is combined with PCA-Net described below.
Kernel methods. A non-parametric technique where the RKHS method is used for encoder, pro- cessor and decoder (Batlle al. 2024).
PCA-Net. A non-intrusive technique with classical POD used as both encoder and decoder, and feedforward network served as processor (Hesthaven & Ubbiali 2018). (Bhattacharya et al. 2021)
POD and oracle. Two POD-based baselines. POD is an intrusive variant of global POD (Volwwein 2013 Oracle is intrusive variant of local POD computed with optimal subspace. In problems consider, error achieved by oracle is the best possible for a given number of basis vectors.
'Recall, that when coefficients in the reduced basis expansion predicted by some model have non- intrusive POD. When basis used generate reduced ODE that later integrated we have intrusive POD.
8

---

Published as a conference paper at ICLR 2026
100
32r
test train
test
oracle
0
train
subspace
S
relative error, %
relative error, %
DeepPOD
DP
POD
P
10
1.5
FFNO
F
DeepONet
D
FFNO,
F
DeepONet,
PCANet
PN
1.5
0.1
kernel
K
DP
Do
DPD: o Fb D PN
(a) Elliptic.
(b) Burgers.
Figure 2: Relative emrors for selected baselines. Label "subspace" refers to subspace regression. For the elliptic problem (a) subspace dimension of ROM methods is bounded by 100, and for DeepPOD, and subspace regression - by 40. Oracle is omitted for the elliptic problem because it has perfect accuracy with 10 basis functions. For Burgers equation subspace dimensions for all methods S 50. FFNO, and DeepONet, refer to an intrusive ROM with bases extracted from FFNO and DeepONet.
Additional results are available Appendix G Here we highlight the most important findings.
Subspace regression leads to competitive accuracy. In Figure 2a and Figure 2b we observe that subspace regression performs similar better than DeepPOD. Among other intrusive methods only classical POD leads to comparable accuracy. Bases extracted from FFNO and DeepONet are worse than global POD in all experiments. Pure regression approaches FFNO, DeepONet, PCANet appear to be less accurate. PCANet similar to kernel methods shows significant overfitting our problems, likely resulting from poor inductive bias of the architecture. Note however, that regression approaches are not directly comparable with intrusive techniques, since they do not require solution of reduced model.
Representations learned by neural networks are highly non-optimal. Representation of all neu- ral networks are inefficient if one compares them with the oracle. For example, by construction of an elliptic dataset, subspace of dimension 10 is enough for perfect accuracy. Both DeepPOD and subspace regression reach eтror about < 1.5% with subspaces dimension 40, DeepONet needs to have > 200 basis functions reach comparable accuracy, and FFNO with 64 basis functions in the last hidden layer lead to 10% relative error. Basis functions built from FFNO and DeepONet are similarly inefficient. The same observations hold for Burger's equation.
4.3 ITERATIVE METHODS FOR LINEAR SYSTEMS
We illustrate subspace regression for iterative methods using D = 2 elliptic problems Said iter- ative methods are deflated CG and two-grid correction for the Jacobi method introduced in Section and explained in more detail in Appendix H Figure 3a and Figure 33b shows average convergence curves on test set and Appendix H contains the rest of relevant data.
Iterative methods are less sensitive to subspace quality. On the training stage, neural networks were presented with data only on first 10 eigenvectors. Despite that, neural networks trained with subspace embedding nearly match the performance of deflated CG with exact eigenspaces of larger size, for coarse-grid corrected Jacobi method convergence speed with learned subspaces is even slightly better. One possible explanation hinted by Theorem 3 is that from distribution of subspaces some information about nearby vectors can be recovered.
Seemingly minor variations in problem setting can lead to substantial variations in the com- plexity of the learning problem. Initially for the Jacobi method we posed subspace regression problem as an approximation of leading eigenspaces of error propagation matrix I-D-1 A, where D diagonal of A. Neural networks with and without subspace embedding completely failed to leam. After inspection the dataset we found that the leading eigenspace contains a complicated mixture of functions with low and high frequencies. Since the learning problem appeared to be completely

---

Published as a conference paper at ICLR 2026
10"
10c
10-1
10-1
relative error
relative error
10-2
10-2
Jacobi
10
M 10
10
M = 10
M - 20
M = 20
M 30
M 30
40
10-41
10-4
50
100
150
200
50
100
150
200
Iteration
iteration
(a) Deflation
(b) Two-grid.
Figure 3: Convergence results for iterative methods. Learned methods are marked with solid lines, and dashed lines correspond to iterative methods with optimal deflation and coarse-grid spaces, M refers to subspace size.
hopeless, we reformulated subspace regression using error propagation matrix damped Jacobi iteration I wD A with 0.9. In contrast to the standard Jacobi method, the damped version leads to subspaces formed by low frequency functions. As evident from Figure [3b] the resulting mapping easily learnable. more detailed report can be found in Appendix H
CONCLUSION
Subspace regression- prediction of subspace from available data is interesting problem with variety of applications including reduced order modelling for partial differential equations, approxi- mation of eigenspaces for eigenproblems, construction of iterative methods for linear problems and optimal control. We formalise subspace regression as a statistical learning problem and introduce several loss functions that suitable for subspace data. For most of the applications considered we observe that the learning problem too complicated even when specialised loss function is used. To simplify learning we propose to approximate given subspace with a subspace of larger dimension. The resulting technique, called subspace embedding, significantly improves accuracy and generalisation gap. The idea of subspace embedding that redundancy typically simplifies the learning process and leads to more robust performance, Even though this strategy clearly helps, introduces large gap between dimensions of optimal and learned subspaces. The same gap observed for the classical operator learning problems, when the neural network is trained to approx- imate solution mapping for parametric PDE. this case learned basis can extracted from the last hidden layer. This neural basis is far from optimal, requiring excessive number of basis vectors to be used for reaching comparable accuracy. Whether this inefficiency in representation can be resolved, remains open problem.
REPRODUCIBILITY STATEMENT
Code used for training, evaluation, and dataset generation is available httts:// com/VLSF/ WVLS//subreg the current version detailed description of architectures, hyperparameters, dataset generation and training details appear in Appendix Appendix [H] Appendix G Appendix
REFERENCES
David Amsallem. Interpolation on manifolds of CFD-based fluid and finite element based structural reduced-order models for on-line aeroelastic predictions. Stanford University, 2010.
10

---

Published as a conference paper at ICLR 2026
Haim Avron, Petar Maymounkov, and Sivan Toledo. Blendenpik: Supercharging lapack's least- squares solver. SIAM Journal on Scientific Computing. 32(3):1217-1236, 2010.
Zhaojun Bai, Patrick M Dewilde, and Roland W Freund. Reduced-order modeling. Handbook of numerical analysis, (3:825-895, 2005.
Pau Batlle, Matthieu Darey, Bamdad Hosseini, and Houman Owhadi. Kernel methods are competi- tive for operator learning. Journal of Computational Physics, 496:1 12549, 2024.
Thomas Bendokat, Ralf Zimmermann, and P-A Absil. A grassmann manifold handbook: Basic geometry and computational aspects. Advances in Computational Mathematics, 50(1);6, 2024.
Kaushik Bhattacharya, Bamdad Hosseini, Nikola Kovachki, and Andrew M Stuart. Model reduc- tion and neural networks for parametric pdes. The SMAI journal of computational mathematics, 7:121-157, 2021.
Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning, vol- ume 4. Springer, 2006.
Ake Björck and Gene H Golub. Numerical methods for computing angles between linear subspaces. Mathematics of computation, 27(112):559-594, 1973.
Barry K Carpenter, Gregory S Ezra, Stavros C Farantos, Zeb C Kramer, and Stephen Wiggins. Dynamics on the double morse potential: a paradigm for roaming reactions with no saddle points. Regular and Chaotic Dynamics, 23(1 2301):60-77, 2018.
Xiangning Chen, Chen Liang, Da Huang, Esteban Real, Kaiyuan Wang, Hieu Pham, Xuanyi Dong, Thang Luong, Cho-Jui Hsieh, Yifeng Lu, al. Symbolic discovery of optimization algorithms. Advances in neural information processing systems, 36:49205-49233, 2023.
Gabriele Ciaramella, Martin J Gander, and Tommaso Vanzan, gentle introduction to interpolation on the grassmann manifold. 2025.
Nicola Rares Franco, Andrea Manzoni, Paolo Zunino, and Jan S Hesthaven. Deep orthogonal decomposition: a continuously adaptive neural network approach to model order reduction of parametrized partial differential equations.
Nicola Rares Franco, Andrea Manzoni, Paolo Zunino, and Jan S Hesthaven. Deep orthogonal decomposition: a continuously adaptive data-driven approach to model order reduction. arXiv preprint arXhv:2404.1884/ 2024.
James E Garrison and llse CF Ipsen. A randomized preconditioned cholesky algorithm. arXiv preprint arXiv:2406.1 11751, 2024.
Magnus Hestenes, Eduard Stiefel, et al. Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards, 49(6):499-336, 1952.
Jan S Hesthaven and Stefano Ubbiali. Non-intrusive reduced order modeling of nonlinear problems using neural networks. Journal of Computational Physics, 363:55-78, 2018.
Jan S Hesthaven, Cecilia Pagliantini, and Gianluigi Rozza. Reduced basis methods for time- dependent problems. Acta Numerica, 31:265-344 2022.
Donald E Kirk. Optimal control theory: introduction. Courier Corporation, 2004.
Andrew V Knyazev. Toward the optimal preconditioned eigensolver: Locally optimal block pre- conditioned conjugate gradient method. SIAM journal on scientific computing, 23(2):517-511, 2001.
Isaac E Lagaris, Aristidis Likas, and Dimitrios I Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE transactions on neural networks, 95):987-11000, 1998.
11

---

Published as a conference paper at ICLR 2026
Robert J Le Roy, Yiye Huang. and Calvin Jary. An accurate analytic potential function for ground- state n2 from a direct- prreettpotentia-ee analysis of spectroscopic data. The Journal of chemical physics, 125(16), 2006.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.0889) 2020.
Lu Lu, Pengzhan Jin, and George Em Karniadakis. Deeponet: Learning nonlinear operators for iden- tifying differential equations based the universal approximation theorem of operators, arXiv preprint arXiv:1910.031 2019.
Brek Meuris, Saad Qadeer, and Panos Stinis. Machine-leaming custom-made basis functions for partial differential equations. arXiv preprint arXiv:2111.0530) 2021.
Brek Meuris, Saad Qadeer, and Panos Stinis. Machine-kkerning-based spectral methods for partial differential equations. Scientific Reports, 13(1):1739, 2023.
Bruce Moore. Principal component analysis in linear systems:
Controllability, observability, and model reduction. IEEE transactions on automatic control, 26(1):17-32, 2003.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686-707, 2019.
Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2003.
Yousef Saad. Numerical methods for large eigenvalue problems: revised edition. SIAM, 2011.
Yousef Saad, Manshung Yeung, Jocelyne Erhel, and Frédéric Guyomare' 'h A deflated version of the conjugate gradient algorithm. SIAM Journal on Scientific Computing, 21(5): 1909:1999, 2000.
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. Advances in neural information processing systems, 33:7462-7473, 2020.
Alasdair Tran, Alexander Mathews, Lexing Xie, and Cheng Soon Ong. operators. arXiv preprint urXiv:2111.138022 2021.
Factorized fourier neural
Lloyd N Trefethen and David Bau. Numerical linear algebra. SIAM, 2022.
UIrich Trottenberg, Cornelius W Oosterlee, and Anton Schuller. Multigrid methods. Academic press, 2001.
Stefan Volkwein. Proper orthogonal decomposition: Theory and reduced-order modelling. Lecture Notes, University of Konstanz, 4(4):1-29, 2013.
David P Woodruff et al. Sketching as a tool for numerical linear algebra. Foundations and Trends@ in Theoretical Computer Science, 10(1-2):1-157, 2014.
Zhi- hii-Qin John Xu, Yaoyu Zhang, Tao Luo, Yanyang Xiao, and Zheng Ma. Frequency principle: Fourier analysis sheds light on deep neural networks. arXiv preprint arXiv: 0et.::901...... 2019.
Yusaku Yamamoto, Yuji Nakatsukasa, Yuka Yanagisawa, and Takeshi Fukaya. Roundoff error anal- ysis of the choleskyqr2 algorithm. Electron. Trans. Numer. Anal, 44(01):306 2015.
Ralf Zimmermann. Manifold interpolation and model reduction. arXiv preprint arXiv:1902.06530 2019.
12

---

Published as a conference paper at ICLR 2026
A PROOF OF THEOREM
1. To show that L1 (A, [3) does not depend on the chosen representative we observe that
k - 2
[A,
= = 11
PB - PAll
(8)
where PA A AT, P BT are orthogonal projectors the columns spaces of and B. When QR decompositions A = QARA. B = QBRp are available, projectors become PA QxQ, PB Q=Qa and identity (8 be verified by algebraic manipulations
2 PP - PAllE - K-P
2 2"h (1QnQ -QQQ)) -Q=QB-QAQ2)) QaQa)) -ー = 1 =2 -(((O) + 22 (r((Q
2
1
11 k -
2 -
(9)
From the equivalent form of loss L (A, B) given in equation one can immediately con- clude that L1 (A, B) does not depend on the representatives A, B chosen from [A]. BB. The reason is projectors are invariant under right GL transformations. For example, PA Px. where A AG and G is arbitrary non-degenerate matrix G E RR
Ā A
A = AG (G A A
AG) G AGG ( (A" G'A' AT. (10)
Now, when we know that (A, B) does not depend on the chosen representatives, it is easy to show that the minimal value of loss is 0 and To see this, select representatives such that QA = (QB 11). where Q3 is block matrix formed from subset of columns of QB and columns of @ are all orthogonal to QB. This selection is always possible since (T - Q=QT) + Q = Representatives selected in this form give
reached when S(B) C S(A).
L, =
= 220.
(11)
where QB E R=xq, q < The last identity follows by construction: QB is composed from columns of QB. Loss becomes zero only if = g. or, equivalently, S(B) S(A).
2. We first show that
(2(A,B;)) min Au Qazl = 1(1 PA) Qnal
(12)
where PA = A (AT -1 AT orthogonal projector on the columns space of A. Using
I Pa) + PA. and (I PA) = (I PA) A we obtain
min ||Au QB²12 = min | - PAQB= (I P^) 2
= min Au - P1 + (I(T PA) 2 = II(T QA) (13)
The last equality holds since PAQB and A share the same columns space. Given that PA does not depend on representative A from [A], and that L2(A, 2) depends on only via PA. we conclude that the same is true for (2(A, B)
3. From equation
I2 we find
E_ [L2(A, B; = E, [[ PA) Q=) E. (Q](11
EE,
E. []= = =I4(A, B). (14)
13

---

Published as a conference paper at ICLR 2026
B PROOF OF THEOREM 2
We provide comments before proceeding with the proof.
In most parts of the text we assumed working with the non-compact Stiefel manifold and in this theorem have data the compact Stiefel manifold (see (Amsallem 2010) for definitions). We specify how compute Q4 and QA having A and A. One may start from any stable version of Cholesky QR, e.g., Gaarrison & Ipsen 2024). Yamamoto et al. 2015), and obtain
Q^ AR=
(15)
where R is Cholesky factorization of Gram matrix A A. i... A R triangular square invertible matrix. To find the derivative of QA need to know the derivative d R. Derivative R can be computed solution to Lyapunov equation
R where a lower
B A,
(16)
after that de R- can be found from Jacobi identity de R-1 R R
In Theorem 2 we use }
- (L) (I)
TIP
r-k 2 to measure the quality of approxima-
tion. follows from the proof in Appendix A that
-- k 2
-
(w) V()),
(17)
where (t),,,, V(t) is loss function defined in Theorem
We can alternatively rewrite
k [(t))(() - (ttV 2
sin (0.)
(18)
using the definition of principle angles 0, between column spaces matrices W() and (tt) (Björck & Golub 1973). Given the later form, it is clear that small values of r correspond better aligned subspaces.
([())(()) (t)
To demonstrate the main result of Theorem 2 we first prove supplementary statement.
Lemma 1. Let A(1) be geodesic Gr(ki,,)) n). A(1) A(I) = Iny One can always construct geodesic B(t), B(I) B(t) Taz on Gr(ki, k2 > bn such that ||000 s | 2 where
[((D) #0.
inequality is strict unless
Proof. Since A(t) is geodesic U sinnnI (T is singular value decomposition of Á(0). Using orthogonality of Á(0) and A(0)
write A(t) A(0)Y cos(2t)Y + U sinn(E))) where
find A(t) = tr()). Without loss of generality we assume that Σn1 * O. Consider
1)
B(t) A(0) ui
(ULUi + 4 cos(o,t))117
9) + Σ sin(o,t)u,
8/
(0)
|
= A(0YY cos Ett T + Us
Et YT
21
(19)
where are columns of Y. u, are columns of U, are diagonal elements of Σ, Σ Σ (01 1)eteT that is, Σ can be obtained from E by replacing En σ by 1. Clearly B(t) is geodesic
and [(() 2 = AA()] -
2
01 < <|AA00] Next we show that principal angles between A(t) and
14

---

Published as a conference paper at ICLR 2026
B(t) all zero. To see this we observe that
(((co0 (x(c0s (2it) (20)) cos [tt + sin (2it) sin (it) (21))) Y
A (t)B(t)=
sin(
= (I- (1 =(((-(--a((ot)aat]) cos(o1))) Y yn sin(ort)) (20)
From the identity above Frobenius norm reads
ln Σ (0,) (0.) (ost((n +
hn
A (AT())(O))
Σ you. + vivi sin'(ont) h1, (21)
and conclude that 0g 0 for all i ,,....,.. ki.
Lemma I also implies that for two such geodesics ||((t)A(t) ||
B(t) b 0.
Now we are ready to show the main result of Theorem 2 split interval of interest E [0.
on subintervals [t,, 1411 of length At. On each subinterval consider three curves: (i) original
continuously differentiable curve V(t) E Gr(k, n.), (ii) approximation of V(t) by geodesic Z(t) E
Gr(k, n.) passing through V(t,) with derivative (tt,), (iii) embedding of ((t) by geodesic W(t)
on Gr(r, n),r > selected as explained in Lemma We start by showing that principle angles between (t) and (t) can be made arbitrary small
2 ||(())(() (W r-k
2
[V()(()) -())Z(11 2(t)Z(t) WOW r-k
(22)
< ((OV()) Z(1)((1T - + (Z0Z))--OW(() 2 (V(()()) Z(1)((1) 1
r-k
2
Now we know that on each interval the distance between (t) and W(t) is bounded by the distance from (t) to the geodesies that passes through (tt) with speed (t,)). Since interval is assumed be small, we expand geodesic Z(t in Taylor series keeping terms proportional (At) and At
and for (I) we use Lagrange reminder (())= = V(.) + V() - 4.), E (4, (4,1, [4,,4]:
[V(()(()) Z(1)(() 1 2(1 4)
(23)
By assumption (t) is continuously differentiable, meaning the expression above can be made bitrary small by selecting sufficiently small intervals [t,, f++1].
To show that derivative of Wtt) can be made smaller than V(t) observe that (t))t z(tt) ||
E on each subinterval where | * 0. Since | 2 [ and ( is
continuous function, we, again, can select sufficiently small intervals such that deviation of |/00
from ((t)) ? on each interval is small enough for | < <|/vcesl? to hold.
C SUBSPACE EMBEDDING EXAMPLE
The proof of Theorem constructive, meaning compute W(t) given V(t) and (t() or its estimation. We select
(S))()(1)) sin(c((t))
(tt) = mt())
ccos(((t) sin( ,((t) 7 cos(2 o(t) x/2 + 7/4ccs(22x) (24)
cos()(tt)
15

---

Published as a conference paper at ICLR 2026
1.0
Gr(3. 1)
_Gr332
0.8
0.6
0.2
(a) 7i (t) E Gr(3, 1)
(b) 1 (t Gr(3,
0.1
0.3
0.4
0.5
(c) Norms of (t) and (2(()
Figure 4: Example of subspace embedding detailed in Appendix C
Curve (l) Gr(1,3) is illustrated in Figure 4a
We next estimate derivatives by splitting on a set of subintervals and taking logarithm on each interval. This derivative is used as explained in Pw - | Y l (7 Pw) we plot (t ((t) in Figure Curve (12(t) appears to be discontinuous, but actually it is continuous owning to Z2 symmetry of compact Stiefel manifold St(1, 3). Norms of derivative are compared in Figure 4c curve 72 manifestly smoother than 7/1.
Lemma 1 to define (t) E Gr(2,3). 3). Since
D PROOF OF THEOREM
In the proof will write F ad G in place of Fj and G2 assuming that & is fixed and the value of evident from the context.
D.I PARTS 1. AND 2.
We order eigenvectors in the increase of eigenvalue E(i == Ai....in Σソ- a which will also call energy this section. To understand how eigenvectors and subspaces are selected for different coefficients Q1..... ap we introduce continuous relaxation of energy E(
- 2,2}. where 2j IR.. In continuous form, surfaces with constant en- ergies are (hyper)ellipsoids of dimension D - 1, so the process of selecting k-th eigenvector or constructing subspace of dimension i can be understood through the following informal algorithm:
1. Select ai...., ap and c 0.
2. Gradually increase c and track ellipsoid Σ/= Σド=0ソ2 =c
3. While increasing add each standard positive lattice point (point with positive integer ordinates) that fall inside the ellipsoids.
4. The order at which lattice points cross an inflating ellipsoid define which eigenvector ap- pears position & and which vectors form eigenspace of dimension k.
To illustrate this process, consider E(21..... 2D)
a + 2222, where we follow procedure outlined above we will see that first lattice points encountered are (1, (2, 1), (3, span of Φ1,1, 62,1, Φ3,1, and the eigenvector that appears on position is Φ3- describe the map from 0i...., aD to du or Vg, this procedure needs to be repeated for all possible positive values of real coefficients ai....,ap.
If
So for considered 01 d2 the subspace of first 3 eigenvectors is
From the algorithm above one deduce that for given aI,... . ap the first time eigenvector with in- dices i1...., iD appears in the sequence of eigenvectors is the first time ellipsoid crosses é1 i...... ip- The position & of this eigenvector will be proportional to the normalised volume of the ellipsoid V.(ai..... a V., where V. V. ap) is a volume D dimensional ellipsoid with semi-axes a1.-..,0 ap ap and V. is volume of D dimensional sphere with radius 1.
16

---

Published as a conference paper at ICLR 2026
(a)
(b)
(c)
Figure 5: (a) Parametric family of ellipsoids passing through point (4,5). (b) Ellipsoid minimal
volume passing through (4, 5). Note that the number of standard lattice points inside is approxi- mately 4
in this case) is asymptotically small for ellipsoids of large volume. (c) Example of non-minimal ellipsoid passing through (4, 5). In the non-minin minimaalll case, the number of standard lattice points inside an ellipsoid passing through a given point can be made arbitrarily large.
20, the error approximation
The first immediate consequence is that d1,..... is piecewise constant function. Indeed, it is clear dI,..., ap can always be perturbed with no change filling order and the change of < < 1, so the eigenvector position & does not change. That proves the first part of the first statement. Next, we need to show that the set of all possible eigenvectors on position is finite.
To see that, we answer the following question: what is the minimal number lattice points one ought to cover with an ellipsoid to reach given lattice point 81....,22 ip. For example, point (1, 1) always reached first. On the other hand, point (2, 1) can be reached arbitrarily late, because one may consider ellipsoids with arbitrary large semi-axis along the second dimension. The position of the point ..... in is known to be the ratio of volumes, so we need to find an ellipsoid with minimal volume that passes through ii, A parametric family of ellipsoids in question and its volume are
&
(ER af
Σ
D
af
1, Ve
π
(25)
[ +1)
D E
(II--1 0z)
See Figure 5a for example of a parametric family in D = 2 passing through lattice point (4, 5). In the expression above used ot to remove positivity constraints. To find minimal volume take derivative with respect to O*
2
OV. 0ak
= 0 = Da
Σ
a
(26)
To find of need compute the nullspace of the linear operator above. It is easy that the solution is af = से. for arbitrary a E R. The example minimal ellipsoid appears Figure 5b
The volume does nod depend multiplicative constant take and obtain normalised minimal volume
min Va (61 -.. iD)
II
(27)
V.
From the considerations above we can conclude that: (i) eigenvector i, i,,...,ip in can not appear on excluding 1,..., appears position k. it can also appear any position > k. Statement (ii) is correct because parametric family of ellipsoids passing through . in contain ellipsoids of arbitrary large volumes unless iz = 1 for all = 1...., D. The example of non-minimal ellipsoid is in Figure
position i unless IT iy k. (ii) if eigenvector
17

---

Published as a conference paper at ICLR 2026
Statement (ii) directly leads to point 2. of Theorem Indeed, since any eigenvector appeared on position & can reappear on arbitrary position / > k. the number of unique vectors that can form low-energy subspace of dimension is the number of eigenvectors position plus eigenvector 1,...,1. That finished the proof of point 1. and 2. of Theorem
Note, that the validity of most of the statements in this section is based assumptions that use continuous relaxation on the problem. particular, we assumed that the position of the eigen- vector proportional the volume of the ellipsoid. Of course in completely discrete formulation this is not strictly the case, but since all statements on the number of eigenvectors and eigenspaces are asymptotic, they will remain valid.
D.2 PART 3.
From the previous section we know the minimal position eigenvector i1 in can appear at. Be- sides that we know that once ....... iD unlocked, it can appear on all positions > k. Give that, the number of eigenvectors on position k reads
D i
#r D)
Σ
Ind
II
S
(28)
ip=l
where Ind is an indicator function. We are interested in asymptotic expansion for large ls and fixed D, the sums above can be approximated by the Euler-Maclaurin formula.
To find asymptotic expansion we will derive recurrence relations for # rk, D). We start by intro- ducing slightly modified function
20 Σ 11=1
D II j=1
# D)=
Ind
i
Sa
(29)
ip=1
Clearly #F (k; D) = k## (k: D) so if we know how to compute # D), we can recover #r(k, D). For D with the help of Bucla-Maclaurin formula we obtain
(, Σ Ind ii S a Σ
Σ 11=1 11
(J,
dx
a a+1
log a+ 1
1
2
2
20
(30)
Next we find recurrence relation
DD413
# D D+
Σ
Σ ip=1
1D41
iD
ip-1=1
ip=1
ip-1
D
(1.
+ (a. D)) (31)
dr
ip41 iD+1
2
is not hard to show that, starting from = 2, recurrence relation can only produce three type of terms: log log)(a), constant term 1. This be seen follows
/
(a/x)
1
1 p+1 log" 1og (a) + (a),
1
la)
2 log
a
dr
c
+
c
+1)
= =clog(k)
+ +
c
c
(32)
2
2a 2
1
a
dr
1.
1
Given that, starting from # r(x, 2) and applying recurrence relations D - times we obtain leading
term
#)
(D- - 1)! log(a) # 1 D) (D - 1)! log(k) P-1
(33)
where last identity follows from the definition of #r D))
18

---

Published as a conference paper at ICLR 2026
D.3
PART 4.
We cannot apply the Euler-Maclaurin formula when is fixed and D is large. To count states under specified conditions will use factorisation prime numbers. For positive integer we can write
p Q1 (p))).().............. 9m, (p) (0),
(34)
where q: (p) are prime factors and o, (p) are their multiplicities. Given this factorisation we can find the number of ways positive integer can be represented as products of D positive integers. All
products of integers correspond to some rearrangement of products in the factorisation of prime factors. The number of such rearrangements is
my
(a, (p) + D 1)!
D), D) = II (D - a).
(35)
This expression is easy to understand if one considers forming the product of D numbers by dis- tributing g-(p) to selected a, (p) among factors for each prime factor gr (r(p), 1,... Mp.
From the expression above, the number of states on position & reads
#(k, D) - Σ (p, D).
(36)
1=1
D is large
mp
(p) +D 1)!
mp
Da-(s)
D (r(
Da(p)
D) II
II
(37)
(D 1)-a)aa(p))
(p))
I1 a, (p))
11 a, (p)
where (2(p) is the prime (big) omega function.
Leading asymptotic expansion of the is the fastest growing term
# (k: D) ma a(p)
Σ
(38)
mi ar (1))
lEarg S a(p)
In the main body of the text we provide a simplified upper bound of this asymptotic expansion. be derived using two upper bounds. First, prime omega function be bounded from above
41 (p) ... 4m (p) (p) > 2 (p) (p) 20(p) 2(p) logz(p).
(39)
Next, the remaining sum can be bounded from above
a (T)! 5 k.
(40)
IT (T)
!Earg max ((p)
These two upper bound combined gives us
#(k; D) < (8),
(41)
D.4 PART 5.
We were unable to compute exact asymptotic expansions for the number of subspaces, our strat- egy in this and next section will be to derive sufficiently strong lower bound by counting selected ways subspaces can be formed.
Consider 2 and = 4. Since surfaces with constant energies are ellipsoids, we can select a large a2 and by gradual decrease of a2 will observe three distinct subspaces:
{Ø1,1, 21,2: $1,3, . (0 21,1 61,2, 21.27 - {Φ2,1 Φ2,2: 01,1, Φ1,2)-
(42)
Similarly, starting from 02 and large decrease of lead to the sequence of subspaces
(d1,1, $2,1 da,1 04.1) (01,2,01, 22,1,, Φ2.13 { $2,2, (21)
(43)
19

---

Published as a conference paper at ICLR 2026
(a) Nu 10
(b) Nab 20
(c) No 40
Figure 6: Illustration of simple greedy subspace embedding technique for elliptic eigenproblem with constant coefficients, See Appendix for details.
Sequences above can be understood as systematic fillings of lattice points along second and first dimensions. First sequence corresponds to filling (0, 4), (1,3), (2, points, and second to (4, 0), (3, 1), (2. In D = the number of distinet subspaces Vi constructed this way equals the number of unordered pairs k1, k2 0 such that k
k.
In arbitrary similarly constructed states can be counted as the number unordered tuples (ki...., kD) of non-negative integers such that E kg k. This is standard counting prob- lem with the answer
(k: + D - -
KD (D- 1)!
(44)
(D 1)K
This provides lower bound asymptotic expansion because our way select subspaces is exhaustive.
D.5 PART 6.
For fixed Is and large D consider eigenvectors with indices (1, (1,..., 1), (1, 11,2,1, 1, 1), (1,..., 1, 2, 1). It is clear that the first eigenvector has indices (I,..., 1) and the rest of them can appear in arbitrary order. This gives at least D(D-1) (D--)--(DD--t. [k- D subspaces.
1)1
E MONTE CARLO EXPERIMENTS FOR EIGENPROBLEM WITH CONSTANT COEFFICIENTS
To illustrate consequences of Theorem 3 we perform simple Monte Carlo experiment. For the case D = 3, we generate 41, 42; 03 from uniform distribution on [0. 1] repeatedly and record distinct subspaces of dimension & = 10. We select unique colour for each subspace and draw them for each
1) This illustration appears in Figure 6a Theorem 3 suggest that the large number of distinct subspaces is a result of selection of k vectors among a small number of candidates. This suggests we may decrease complexity of the function from coefficients to subspaces by predicting excessive number of vectors. In our experiments with neural networks this redundant mapping learned, here we build the mapping with simple greedy strategy. In place of function a @22 d3 - Vio we consider 01, 02, Vio {VI,..., vK where U1...., Ey are first & most abundant eigenvectors. By appending additional vectors we decrease the number of distinct subspaces. The result of this greedy simplification appear in Figure 6b with 10 additional vectors and in Figure (6c with 20 additional vectors. As we see the number distinct subspaces rapidly decreasing.
point 01, 02 ag on the plane perpendicular (1
20

---

Published as a conference paper at ICLR 2026
F DETAILS ON NUMERICAL EXPERIMENTS FOR EIGENPROBLEMS
F.1 DATASETS
We generated two datasets for D elliptic eigenproblem with uniform Dirichlet boundary condi- tions
div k grad Φ; A l|obil2 1, x E [0,1]²
(45)
For both datasets we used uniform grid 100 100 and unite-difference discretisation. Components of diffusion coefficients were generated from the same distribution for both datasets. Diffusion coefficient is generated as follows:
Gaussian random field ub is generated from N (0, (id (id-124) (()-)>))),>> 7 1 r- 2- 20
2. Diffusion coefficient is computed
+ (3 ax) (tanh (sab) + /2 with 1,
B = 50, 8 =
For one D 2 dataset h hi and for another k # 12 but both are i.i.d. random fields generated as described above. In the main text only results for h = 22 are reported.
For D 3 elliptic eigenproblem use setup analogous D 2 but grid of size 30 30 30 and Kı = 22 33 generated the same way as explained above with parameters Too 50, 1.
r
2
For QM problems datasets are defined by distributions for potential functions.
For D we use
( 1 1- cxp
((1))))
Q1
<Te
V(r)
Te To
p(r) +1
(46)
Q2
r 2
where
(x) x x
(z) -1---1
-
(47)
x+1
and q (x) is polynomial of degree deg. Polynomial (2(x) has the same form.
In dataset is by selecting uniform grid with 100 points on the interval [0. 10), re is sampled from uniform distribution on the interval [1, 8). d sampled from uniform distribution on the interval [10, 40), both 9I and O2 has order 10, for qı all coefficients (including c) are sampled from uniform distribution on [0, 5), for ø2 coefficients of are sampled from uniform distribution on [0, 10] and sampled from uniform distribution on the interval [1, 11].
For D 2
V(z,8) ()=,))=VV V
V (1 cu)² + (u) + V/2 V(x+ ou)² (y+ cc)3
(48)
where are component of random normalised vector, V2r.. Potentials V, and 1/2 are i.i.d. with parameters: order of polynomial is 2, re uniformly distributed on [1, 5). d is uniformly distributed on [10, 40), coefficients of q are sampled from uniform distribution on [0,3] and is sampled from uniform distribution on [10, 13). To discretise the problem we use finite difference and uniform 100 x 100 grid on [-7 7,7 7|2
F.2 ARCHITECTURES AND TRAINING
In all cases we used FENO architecture, with GELU activation functions, that is completely specified by: number of layers N Nayars numbers of features in hidden layer Nisauues number of Fourier modes in spectral convolution Nmodes. Since all our loss functions are scale-invariant, the output of FFNO architecture was normalised.
We use Lion optimiser (Chen et al. 2023 with weight decay. Parameters of the optimiser are learning rate Ir, rate decay factor Ydecay and number of transition steps Namaay-
For D 2, eigenvalue and QM problems, and also D = QM problem we perform grid search with parameters: Nayers [3, 4, 5), N Nieatures 64, Nmudes [10, 14, 16), Ir [10-3, 10-2,, 'deaay
0.5.
21

---

Published as a conference paper at ICLR 2026
Naecay E [100, 200). batch size fixed to 100, number of train samples is 4000, number of test samples is 1000. Architecture is training to approximate subspace spanned by first 10 eigenvectors. Number of epoch is 1000. When architecture is trained to predict individual eigenvectors, the same grid search applies.
Since architectures used neural operators, by first training architecture on coarse grid 32 x 32, and then retraining 3 most accurate architectures (according the results of grid search the coarse grid) of the fine 100 100 grid.
able to speed up grid search by factor of
For = 3 grid search not practical, so we select Nayers = 4, Nieaturs = 128, Nmodes = 16,
= 10 10-3, "ieeay = 0.5, Ndecay = 100. The size of the train set 800, the size of test set is 200. The number of epochs is 1000. Architecture is training to approximate subspace spanned by first eigenvectors.
F.3 ADDITIONAL RESULTS FOR ELLIPTIC EIGENPROBLEMS
Additional results are available in Table 3 Table 4 Table 5 best observed result among three best grid search runs.
Results in brackets indicate worst and
Table 3: Comparison of (A, B) and L2 B) loss functions for kn k2
L1(A B)
L2 B;
Nauh
Eqmin
Evest
4min:
Errain
Evst
Urain; 8
10
00.216, 0.235]
[0.292, 0.314
4124 + 443
[0.244, 0.252]
[0.296, 0.302]
38951:524
20
0.038, 0.048
0.046, 0.052
5962 +
[0.046, 0.052]
[0.049, 0.058)
4074+1215
30
0.024, 0.029
0.028, 0.033
7902 + 85
[0.026, 0.033]
[0.029, 0.037]
3973+
40
0.018, 0.025
0.021,0.029)
10842 570
[0.017, 0.024]
[0.02, 0.027]
4270++95
Table 4: Comparison of L, B) and (2,,, 2; loss functions for K1 & 12
T. (A, B)
Tz A, B;²)
Naub
Eyrain
Evsst
tuain:
Euain
Evst
tuain
10
00.305, 0.377]
[0.407, 0.421]
4035 + 295
[0.312, 0.335]
[0.386, 0.389]
3661 4:018:308
20
0.066, 0.089]
0..09, 0.10
5817 + 219
[0.092, 0.092]
[0.103, 103]
4262+
30
0.055, 0.05]
0.063, 0.063
7966 + 84
[0.042, 0.05]
[0.052, 0.059]
19911168
40
[0.035, 0.036]
00.045, 0.047]
11019 + 252
[0.034, 0.038]
[0.041, 0.046]
4238 112
Table 5: Zz-adjusted L2 loss.
11
K2
1 # Kz
Neig
Firai
Fiest
firain $
Farain
Fiest
fuain
0
[0.009, 0.012
[0.036, 0.038
3666 + 476
(0.014, 0.028
(0.068, 0.07]
3280 + 60
1
00.08, 0.042
(0.158, 0.165
3312 + 77
0..031, 0.033
(0.19,,
0.218
3956 + 162
2
[0.046, 0.048
[0.359, 0.373
4168 + 82
0..45, 0.053
(0.553, 0.563
4140 + 87
3
[0.046, 0.057
00.41,005555
4168 + 81
0.054, 0.066
0..74,, 0.779
4101 + 87
4
00.068, 0.084
0..554, 0.769
4007 + 204
0.002, 0.078
0.945, 0.97]
4101 + 87
5
[0.073, 0.075
00.897, 0.905
4041 264
0.084, 0.094
1.087., 1.098]
4090 - 83
22

---

Published as a conference paper at ICLR 2026
F.4 GENERALISATION DIFFERENT GRID SIZE
For all subspace regression problems we use FFNO. Since FENO is neural operator it should be discretisation agnostic. Here we report results for model trained grid 100 on D quantum mechanics eigenproblem and tested on grids of higher resolution. For each grid we generated new test from the same distribution as specified in Appendix F.1
We see approximately linear increase of relative error with resolution. Owning good initial accuracy on grid N2 100 the relative error remains under 10% for grid with N₂ = 500.
The results are available in Table 6
F.5 SMOOTHNESS OF NEURAL NETWORKS TRAINED WITH SUBSPACE EMBEDDING TECUNIQUE
To empirically measure smoothness of learned map, introduce several "smoothness indicators":
1. Taylor indicator
12;11 W +12) - - V +12) (f2)
(49)
[W(G - 1/2)11
Taylor indicator a relative error of linear model. One expect that: (i) unless i linear function, when increases relative error also increases; (ii) when smoothness increases Taylor indicator decreases.
2. Average cosine
Caf -
Ecoss cos cos tf)
(50)
where cos, (A, B) are cosines of principle angles BBööckk Golub 1973 and is target at point f1. For smoother maps average cosine increases until it reaches maximal value of 1.
3. Frobenius norm of the directional derivative
1 D
d di
Flf1. f2:
N(f++
1122
(51)
1=0
where D the subspace size. The magnitude of the directional derivative computed by automatic differentiation and it is expected to decrease when smoothness increases.
Each indicator depends on two features /1 and f2 and real number E [0, 1]. Results reported for neural network trained to predict eigenspaces for elliptic eigenproblem with k k2 and K1 * k2 (see Appendix F.T for description) with loss func- tion L2(A, B: 2). For each indicator we provide values for several I averaged over 1000 randomly selected feature pairs f1. /2. Results are reported in Table 7 Table 8 Table
Table 7: Frobenius norm of the directional derivative smoothness indicator (51).
N subspacc
K1
k1 + 122
10
10.73
8.26
20
9.72
8.15
30
8.78
8.14
All indicators clearly demonstrate the improve in smoothness when the size of embedding N subspacc increases. Results for average cosine indicator (50)
40
8.26
6.88
Table 6: Network is trained on resolution N₂ 100 for D with increased resofution.
QM problem and evaluated grids
N test error, %
100
150
200
250
300
350
400
450
500
0.83
1.37
2.12
2.95
3.96
4.32
5.33
5.32
6.79
23

---

Published as a conference paper at ICLR 2026
indicate that learned mapping effectively average information about subspaces for distinct features. It becomes especially clear if one compares results for ! = 1.0 with average cosines computed between targets D ش cos; t2): for k1 k2
average cosine is 0.51; for + k2 average cosine is 0.53.
Table 8: Taylor smoothness indicator (49).
k
12
ky # h2
N subsspace
= 0.25
= 0.5
= 0.75
= 1.0
1= 0.25
= 0.5
1= 0.75
1= 1.0
2.41
5.1
7.48
10.71
1.8
3.87
5.99
8.21
20
2.12
4.57
7.05
9.68
1.72
3.77
5.86
8.08
30
1.88
4.09
6.34
8.73
1.71
3.77
5.87
8.08
40
1.76
3.84
5.98
8.22
1.42
3.15
4.94
6.8
Table Average cosine indicator (50
k1
h1 7 12
N Nauhspace
1 = 0.25
[= 0.5
- 0.75
= 1.0
== 0.25
1 = 0.5
1= 0.75
/= = 1.0
10
0.81
0.74
0.67
0.51
0.82
0.74
0.67
0.53
20
0.94
0.88
0.82
0.68
0.94
0.88
0.83
0.72
30
0.97
0.93
0.88
0.78
0.97
0.93
0.89
0.81
40
0.98
0.95
0.91
0.82
0.99
0.96
0.93
0.87
F.6
SUBSPACE REGRESSION COMBINED WITH LOBPCG
Subspace regression be also used to improve results for classical iterative eigensolvers. We demonstrate this for LOBPCG, which is a matrix-free iterative eigensolver that can approximate extremal eigenspaces (Knyazev 2001). A notable feature of LOBPCG is possibility of hot start: when approximation to eigenspace of interest available, can be used at the initialisation to speed up convergence. In Table 10 we report such speed up for elliptic eigenproblems described in Appendix [.I Metrics in Table Th computed for test set, and maximal number iterations for LOBPCG is set to 1000. Convergence plots are available in Flgure 7a and Figure [bb It is evident that initialisation by subspace regression lead to both smaller number of iteration and better final error. Note, that the overall cost of method is dominated by the cost of iterations.
Table 10: Performance of LOBPCG with initialisation by subspace regression compared with ran- dom initialisation, N₁, number of iterations until convergence.
h
12
k1 * l2
initialisation
N Suuhquace
N₂
relative + std
N₂
relative error + std
subspace regression
10
410 245
0.093 1.109
332 + 213
0.028 + 0.179
subspace regression
20
288 + 197
0.159 + 4.12
249 + 174
0.02 + 0.087
subspace regression
274 209
0.047 0.554
227 + 172
0.022 + 0.119
subspace regression random
40
263 + 200
0.131 + 3.218
221 + 160
0.02 + 0.089
685 265
7.86 + 30.16
609 + 231
11.42 + 35.49
24

---

Published as a conference paper at ICLR 2026
10a
subspace init
subspace init random init
random Init
reisdual, L2 norm
109
reisdual, L2 norm
10a
10
10-2
2 2-10-5L
5-10-5
200
400
600
800
0
200
400
600
800
iteration
iteration
(a) kn
k2
(b) D 2, k1 + kz
Figure 7: Convergence plots for LOBPCG with and without initialisation by subspace regression.
G DETAILS ON NUMERICAL EXPERIMENTS FOR PARAMETRIC PDES
G.1 DATASETS
For the D = 2 stationary diffusion equation we reused datasets described in Appendix To generate forcing terms for each we select 10 eigenvectors 1; l,..., corresponding smallest eigenvalues and compute exact solution u Σ, di2: where ~ (O0,,) 1). Forcing term f(x) corresponding to this solution is Σ, 0:24.
For D + Burgers equation
sample random diffusion coefficient v(x) and initial condi- tion un (z). In both cases we first sample Gaussian random field from A (0, (id "x) To sample 40, 4. Diffusion coefficient is computed For initial conditions we use random field with a = 10, 2. Uniform grid is used discretise equation with 128 points for and 64 points for t. Time interval is [0. 10-1] and E [0, 1].
random field v (x) for diffusion coefficient
as V 5 10 + (1 + vanh(30y))
G.2 ARCHITECTURES AND TRAINING DETAILS
We performed grid search for all methods that involve learning. We start by describing hyperparam- eters of architectures.
FFNO used for subspace regression, standard regression, DeepPOD, and for intrusive techniques with basis extraction. Parameters of FFNO used for grid search are described in Appendix F
DeepONet is used for standard regression, and the intrusive technique with bases extracted from branch net. As branch net of DecpONet we select a classical convolution network with downsam- pling by factor of 2 along each dimension and the increase of the number of hidden features by factor of 2. Branch net is defined by the number of features after encoder Ne.b. kernel size of convolution kg, and number of layers Ng. Trunk net is MLP which is defined by the number of hidden neurons Ny.,, number of layers N₂ and the size of basis on the output layer No. In out grid searches we used Ne, E [4, 5), K, E [3, 7J. N₂
N₁.- f. = No E [100, 200), N, E [3, 4].
PCANet is defined by the number of POD basis functions used to compress feature and targets Np and Np.e. Number of MI.P layers NMLP and hidden units N, M.P.P. our experiments we use Np.j E [100, 200, 300, 400, 500] for elliptic equation and Np.f [50, 80, 100] for Burgers equation, No.t E [100, 200, 300, 400] for elliptic equation and Np.t E [100, 150, 400] for Burgers equation, NMLP E [3, 6] for elliptic equation and NMLP E [3, 7] for Burgers equation, Np.MLP
[100, 200, 300, 400, 500] for elliptic and Np,LLP E [100, 300, 500] for Burgers equation.
Hyperparameters of kernel methods are the type of kernel and the number of POD basis functions used to compress features and targets Np,, and Np.-- We use Matern and RBF kernel and Np.f E [50, 100, 150, 200), Np.t e [50, 100, 200)
25

---

Published as a conference paper at ICLR 2026
95
82
POD
POD DeepONet
DeepONet
FFNO
FFNO
relative error, %
DeepPOD
relative error, %
DeepPOD
16
subspace regression
subspace
regression
4
2
10
50
100
150
200
10
50
100
150
200
Nbasis
Nbasis
(a) D = 2, ly = K2
(b) D = 2, k1 * kz
Figure 8: Relative errors for two stationary diffusion equations depending on the number of basis functions for several basis construction methods. Bases constructed with DeepPOD and subspace regression lead to the most accurate intrusive methods.
To train neural network we used Lion optimiser with Ir [5 10-5, 10-1] for FFNO and Ir [10- 10 10044] for all other architectures, "deeay = 0.5, NNecaa E [100, 200]. We use batch size 10, train PCANet for 3000 epoch and other networks for 1000 epoch.
G.3 ADDITIONAL RESULTS
We provide two additional results. For elliptic equations compare optimality of learned bases. The results are in Figure 8a and Figure 8b) Interestingly, global POD) leads to a better basis than both FFNO and NeepONet. For DeepONet this is expected, since trunk does depend on parameters of PDE. For FFNO the result is more surprising, because the basis is extracted from the last hidden layer, so explicitly depends on parameters. Bases of DeepPOD and subspace regres- sion are the most optimal one, still they are highly non-optimal when compared with local POD. For the Burgers equation we compare Deep- POD and two variants of subspace regres- sion in Table 11 SubReg 100) subspace regression trained to approximate subspace spanned by first 10 local POD basis func- tions and SubReg(5) was trained with 5 lo-
Table 11: Relative errors for Burgers equation. Target for (ub is subspace of dimension n.
N Subspace
DeepPOD
SubReg(10)
SubReg (5)
cal POD basis functions. DeepPOD is an unsupervised method and trained with
10
16.37%
22.79%
14.34%
20
10.0%
11.73%
10.68%
the whole trajectory. This implies meth-
30
5.42%
6.29%
5.03%
ods are sorted from left right in the de- crease of information they receive about
40
2.57%
3.46%
3.01%
50
1.8%
2.49%
1.74%
solutions. Interestingly, the SubReg(5) method, learning the smallest subspace performs better almost uniformly. A possible explanation is that an optimal subspace of dimen- sion 5 leads to good enough accuracy and is easier to learn than larger subspaces.
H DETAILS ON NUMERICAL EXPERIMENTS FOR ITERATIVE METHODS
H.1 DEFLATION
Since for a deflation problem one needs to approximate eigenspace spanned by eigenvectors with small eigenvalues, reused dataset and network trained for elliptic eigenproblem. The description of training and datasets is available in Appendix Fl Results the main text are for k = 12 for kn * K2 convergence plots given in Figure
26

---

Published as a conference paper at ICLR 2026
10
10-1
relative error
10
10
CG
10
M
=20
M = 30
10-5
50
100
150
200
iteration
Figure 9: Convergence results for deflated CG, elliptic dataset D 2 with kj o k2. Learned methods are marked with solid lines, and dashed lines correspond to iterative methods with optimal deflation space, M refers to subspace size.
H.2 TWO-GRID METHOD
Elliptic equation dataset We consider 2D elliptic equation on the unit square Ω = 1)2
with homogeneous Dirichlet boundary conditions (5). We aim learn mapping k(x) + S(V). Variability of the dataset comes from the spatially heterogeneous coefficient function k(x)).
Each sampled coefficient function is strictly positive random field built in three steps:
1. Draw i.i.d. Fourier coefficients on a square index set K = {o.... M - 113 and form a real field by summing complex exponentials. We additionally introduce Fourier-space weight Wkk = (1 + X (k())1 to control the high-frequency components.
2. Multiply the real field from the previous step by A2, then apply a hyperbolic tangent func- tion to control the contrast of the coefficient field values.
3. Rescale the field to the prescribed interval [a, 3] to ensure strict positivity and enforce controlled contrast ratio of 3/a.
Exact procedure to generate the 2D field is:
so v) Re
Σ
1 A1 k ck N00
(O......M111
s(x, y) = tanh (A2 so(x,))).
k(x,v) (B a)) s(x, v) +
k(x, v) E (a,
2
The equation is descritized on uniform grid with point finite-difference stencil, yielding sparse, symmet- ric positive-definite matrix. One can observe sampled normalized coefficient function in Figure 10
0.30
Target subspace We aim to predict a coarse-grid sub- space for the two-grid method, which applies coarse- grid correction
- x V(V -V Ar)
with weighted Jacobi smoothing
8.20
wD-1 As))
Figure 10: Sample coefficient function.
27

---

Published as a conference paper at ICLR 2026
before and after. This coarse projection V is learned troblem-specifte subspace from the coefficient field. In this setup, our projection matrix V spans the leading eigenspace of the error propagation matrix I - wD- 1-D--AA. A. While the relaxation parameter w used to quickly dampen the fast modes, the coarse-grid projection removes the low-frequency components, resulting in rapid overall convergence. Thus, target subspace for regression consists of the first few eigenvectors of / - wD= A, sorted by the absolute values of their eigenvalues.
Subspace relaxation Initially, we experiment with under-relaxation (i.e. w 1.0) in the Jacobi smoother. As expected. this leads to a target subspace that contains both high- and low- frequency modes. Neural networks are known to have a spectral bias towards low frequencies. In FNO-type models, this spectral bias is especially pronounced because these models truncate high- frequency modes in the Fourier domain. We address this by reducing the relaxation parameter to w 0.9, which makes the leading subspace components dominated by slow modes. It is worth noting that this adjustment is consistent with the respective roles of smoothing and coarse-grid correction discussed above. In Figure [ can observe the first five eigenvectors of the target subspace for representative sample under both relaxation settings (w 1.0 and w = 0.9).
Figure 11: First five eigenvectors of the error propagation matrix / wD-1A. Top: = 1.0. Bottom: w 0.9.
Subspace prediction Before applying the neural network's predicted subspace, we perform QR decomposition the predicted matrix W to obtain an orthonormal basis Qw. We assess the quality of the predicted coarse subspace in the two-grid method using three metrics:
1. Cosine angles between the true subspace V and the predicted subspace, computed as the singular values of Q V. A value closer to 1 indicates better subspace alignment.
2. Relative reconstruction error e = min, ||V W ull2 for each true basis vector V, computed
[I(] Qw (Tr V ||2. A smaller value indicates that the predicted subspace reconstructs Vi more accurately.
3. Two-grid convergence rate, measured by the spectral radius p of the two-grid iteration operator T', We estimate via the power method by repeatedly applying T' to a vector:
Tuk
U
T
A smaller spectral radius indicates faster asymptotic convergence.
In Table 12 we report these metrics for the best-performing models and for the ground-truth target subspace. Across all experiments, the predicted coarse subspaces achieve slightly better two-grid convergence (i.e., lower spectral radius) than the ground-truth subspace. These results are partic- ularly interesting since smaller subspaces yield rather high cosines and reconstruction errors. The increase of size of the predicted subspace improves the quality of the reconstruction and does not degrade in effect on iteration compared to true exact subspace. It is also worth noting that both training objectives yield similarly effective subspaces.
28

---

Published as a conference paper at ICLR 2026
Subspace prediction metrics for two-grid method and Jacobi iterations. We report av- eraged values over the test set. Methods' column values: Exact subspace uses the exact leading eigenvectors of eтror propagation matrix (ground truth). loss and L2 loss denote subspaces pre- dicted by F-FNO trained with the respective objectives. The Jacobi baseline reports the spectral radius with no coarse correction (w 1.0).
Table 12:
Subspace size
Method
Cosine
Rec. error
Spectral radius
Jacobi
0.9976
Exact subspace
0.9917
10
L1 loss
0.845
0.411
0.9910
L2 loss
0.859
0.392
0.9908
Exact subspace
0.9858
20
I loss
0.960
0.222
0.9852
L2 loss
0.962
0.217
0.9852
Exact subspace
0.9799
30
Lj loss
0.986
0.140
0.9790
L2 loss
0.986
0.139
0.9790
Exact subspace
0.9745
40
L1 loss
0.994
0.097
0.9730
Lg loss
0.994
0.097
0.9731
Data and training details We generate two different datasets with 32 and 100 interior grid points.
Both datasets use MM = 100 Fourier modes, Å = 0.1, A2 = 1, interval [a, 3] = [1, 50). and w 0.9 in error propagation matrix. Each dataset contains 1, 000 training and 200 test samples. While
neural network predicts subspace of sizes {10, 200,20, 30, 40), target subspace always contains 10 basis functions.
We train the Factorized Fourier Neural Operator (F-FNO) model (Tran et al. 202I). We first conduct extensive hyperparameter search the 32 dataset with:
Number of retained modes: {10, 14, 16).
Number of processor layers: (3, 4, 5).
Learning rate: {10 - 10-4..
Step-decay every {100, 200) epochs.
By default, the batch size is 64, training runs for 1, 000 epochs, and each processor layer has 64 features. We repeat this search for both La and T-2 losses and for predicted subspace sizes (10,20) 30, 40). We then select the top-3 hyperparameter configurations per subspace size and loss (by two-grid spectral radius) and train on the dataset with 100 grid points. Throughout the paper, we report results for the best configuration on the dataset with 100 grid points.
29

---

Published as a conference paper at ICLR 2026
1 SUBSPACE REGRESSION FOR OPTIMAL CONTROL
We consider optimal control of D heat equation with homogeneous Dirichlet boundary conditions
00(x,t) = div k. grad (z) b(z) + Σ wy (z) (t),,(t,
0t
= =(Vi,
(52)
A min L UU)) ((T) + =22// L u(t) u(t))u(t). u(t)
The problem has simple interpretation. With no control for sufficiently large T system reaches steady-state, which can be computed as solution of the linear system div ggaa ((x, t) b(x)). The objective function contains the term (7) (T) 2 O. Optimal control minimises amplitude of (TT). Given that y() is a projection of state variable & on vectors vh, the result of optimal control
reach with E, (d. th)) 2bs where 1y is any basis in subspace spanned by y/s.
Control problem 52 is not linear-ue linear-quadratie regulator in its standard form since b(x:) is present. To get rid of 6(x) and use exact solution for quadar-qu regulator we first discretise PDE (52) and after that introduce additional (constant) variable
((()) = 4 )) dt ((())
((1) ((t))
+ -00)
W
((0) = Φo ((0)
",
((xx))
(53)
min L 1 (T) 2
T Vo ((T) A 2 u(l) u(().
Problem 53
is a standard linear-quadratic regulator, so the exact form of value function is known Optimal control can be computed follows
Kirk 2004
Či Cn AC X CW, C12 0, Ci((T) 0,
Čn + C11 + AC - -- = 0, CT((T) = vv
(54)
uu() = - (C11 C12b(x)
As model reduction method we apply balanced truncation closely following (Moore 2003): solve Lyapunov equations to find controllability and observability Gramians, (ii) find eigendecomposition of controllability Gramian and select coordinates system where controllability Gramian identity matrix, (iii) find eigendecompositton of observability Gramian and select coordinate system where controllability and observability Gramians coincide, (iv) from the composition two coordinate transformations build degrees of freedom corresponding to largest eigenvalues of both controllability and observability Gramians.
1.1 DATASET
We use random gaussian random field A (0, (id 10 aA))) to generate (x), (x)),, diffusion coef- ficient k(x), initial conditions up(x) and forcing b(x)). For wi, 1h,, uz)), b(x) take 5 and n = 4, t and w; are further orthogonalised with QR, 30 i.idd. Wi and thy are generated for each dataset sample; for diffusion coefficient we
6, n and process generated random field similarly to Burgers equation k(x)) = 5 x 10-3 + (1 + tanh(5x))/10. Equation discretised uniform grid 128 128, x E [0, 1). [0,5]. Dataset consists of 1200 samples, 1000 for train, 100 for validation and 100 for test. Optimal reduction by balanced truncation is computed for each sample and later used for subspace regression.
1.2 ARCHITECTURE AND TRAINING DETAILS
We FFNO and precisely the same grid search as for the Burgers equation.
30

---

Published as a conference paper at ICLR 2026
Table 13: Results for control.
N basis
exact
(A B)
(2(A, B3z)
E.
E.
E,
Eo
E.
E.
4.16%
4.07%
12.56%
12.21%
10.91%
10.59%
20
8.8%
8.47%
7.96%
7.65%
30
9.49%
8.94%
8.74%
8.42%
40
7.26%
6.99%
7.41%
7.15%
50
7.2%
6.94%
7.09%
6.88%
1.3 RESULTS
Results are summarised in Table 13 on first 10 basis vectors obtained with balanced truncation. As metrics we use relative observa- tion eтror E,, at time T and relative full state error E, at time One can observe that subspace embedding techniques improve accuracy for both loss functions. Interestingly. La leads to slightly better error for small subspace sizes. Overall accuracy is acceptable but does not reach optimal performance reported in the first columns.
We train neural network with two subspace regression losses
SENSITIVITY TO TRAIN-TEST SPLIT AND THE CHOICE OF HYPERPARAMETERS
All results the main text are reported without error hars, that are usually computed by varying ran- dom initialisation train-test split. All our experiments involve extensive hyperparameter search, so computing statistics for distinct initialisations or train-test split would require order of magnitude more compute. Here demonstrate for elliptic eigenproblem (Appendix F.I) and Burgers equa- tion (Appendix G.I that sensetivity to train-test split is relatively small and not significant to main conclusions reached by analysis of "single-run" results.
Table 14: Sensitivity to train-test split for elliptic eigenproblem and Burgers equation, subspace regression trained with L2 (A, B; 2) loss function.
elliptic, k1 = K2
Burgers
N Naahquace
train error + std, %
test error + std, %
train error + std, %
test error + std, %
10
24.37 + 1.47
30.38 + 1.34
25.5 + 1.17
25.71 + 1.72
20
5.69 + 0.53
6.46 + 0.26
11.15 1.2
11.0 1.43
30
3.34 + 0.26
3.79 + 0.35
7.04 + 0.35
7.01 + 0.64
40
1.84 + 0.12
2.12 + 0.11
4.56 + 0.32
4.72 + 0.61
Table 15: Sensitivity to train-test split for elliptic eigenproblem and Burgers equation, subspace regression trained with L2 A, 2) loss function.
elliptic, k hz
Burgers
train error
test error [min, max], %
train error
test error [min, max], %
N Subspace
min, max], %
[min, max], %
10
23.81,
28.67, 39.05)
222.01, 32.64
[23.79, 34.62)
20
4.74, 100.0]
5.51, 100.0]
9.53, 20.14]
10.06, 21.05]
30
2.77,8.6)
3.33, 9.06)
5.58, 21.81]
5.84,
40
[2.06, 7.05]
7.34
4.3, 43.43
(4.11, 283.33)
31

---

Published as a conference paper at ICLR 2026
For both Burgers and elliptic equations randomly split dataset train and test set 5 times and report mean metrics and standard deviation. Networks are trained for hyperparameters found by grid search. The results available in Table 14 Variability is clearly present, but it is not pronounced enough.
Variability train-test split should be compared to sensitivity to the selection of hyperparameters reported in Table [5] We see that for certain cases the span of best and worse performance reaches 100% which is much higher than variability train-test split, Note that for elliptic equation grid search was performed on downsampled dataset with 32 x 32 grid.
32


Reviews:
[{"review_id": "McY8l2mx4d", "reviewer": "Reviewer_LUwp", "summary": "The paper introduces Subspace Regression, a learning framework in which the target is a linear subspace rather than a vector or a scalar. It arises naturally in reduced-order modeling (ROM), eigenspace approximation, deflation methods, and parametric PDEs. It formulates the subspace regression on the Grassmann manifold $Gr(k, n)$, proposes two new loss functions invariant under right $GL(k, V)$ actions, and proposes the subspace embedding technique, which predicts a higher-dimensional subspace that contains the target, theoretically shown to smooth the mapping and reduce complexity. I mainly focused on checking the proofs of the theorems and the examples, which look very rigorous to me.", "questions": "1. Just a tiny suggestion: Sec. 2.1, I suggest using [W] to denote the equivalence class. \n2. Line 88, why the parameteric model maps to Gr(r, n), r \\geq k, which is more expressive than needed?\n3. In the Schrödinger eigenproblem, I do not really get why it is a subspace regression problem. Do you wish to learn the coefficients or the basis $\\phi_i$? Also, in line 113, what is $f$? $\\alpha_i$ should be complex-valued.\n4. How is the embedding dimension $r > k$ chosen in practice? \n5.  How sensitive is the training to the choice of r.v. in the stochastic loss $L_2$?\n6. For high-dimensional PDEs, does the cost of orthogonalization or QR decomposition outweigh the benefits of manifold losses?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "1. Formally consider the subspace regression problem with differential geometry. \n2. The list of example problems (Sec. 2.2) to solve with subspace regression is very valuable and illustrative, IMHO. \n3. I personally like the writing quite a bit, which is concise and right to the point; not much redundancy in the sentences. \n4. Theorem 1–3 rigorously justify the loss functions and the subspace embedding principle using differential geometry and asymptotic complexity arguments.\n5. Comparisons with classical interpolation, FNO, DeepONet, and POD variants show consistent advantages, especially in high-dimensional parameter spaces.", "weaknesses": "1. Tiny suggestion: for QR factorization, it is more precise to say $A = Q_A R_A$, where $Q_A\\in \\mathbb{R}^{n\\times k}$ since for a tall matrix $A$, you can factorize it into a unitary $Q \\in \\mathbb{R}^{n\\times n}$.\n2. Another small suggestion, in Thm. 1, for $L_1$, if $A = B \\cdot G$, for some $G\\in GL(k, V)$, then $L_1(A, B) = 0$, which explains the $p$ term in the loss, which might look a bit weird at first look. \n3. Yet another thing: for the parametrization of $V$, you wrote $V(r) \\colon \\mathbb{R}^p \\rightarrow Gr(k, n)$ (which I understand), which, however, confused me at the first glance: I thought it mean for each $r$, $V(r)$ is a map from $\\mathbb{R}^p$ to $Gr(k, n)$. Better to write $V \\colon \\mathbb{R}^p \\rightarrow Gr(k, n)$.\n4. Thm. 2 shows that larger subspaces help with reducing the magnitude of the tangent vectors, but it doesn’t quantify why it helps with optimizing the loss function, e.g., I think Lipschitz analysis of the tangent vectors might help. \n5. Many results are reported as single-run metrics without confidence intervals or error bars.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "kbVk0wswGm", "reviewer": "Reviewer_fKUi", "summary": "The paper studies a subspace regression problem that aims to learn a mapping from a parameter set to a manifold of subspaces. One central challenge is that the target of the regression is itself a *set of spaces*, rather than points in a Euclidean domain. The main contributions lie in designing appropriate loss functions and embedding techniques that allow these spaces to be compared, optimized over, and represented effectively within a learning pipeline.", "questions": "There are several technically unclear points in the presentation:\n\n1. Line 209: the meaning of “interpolation” is unclear and should be defined precisely in this context.\n\n2. Theorem 3: it is not explained why the functions $F_k$ and $G_k$ are piecewise constant, nor why they take distinct values. A justification or reference is needed.\n\n3. The description of the **SUBSPACE EMBEDDING** technique is very unclear. How does it work in practice? Does it imply that the learnable data consist of higher-dimensional subspaces, and if so, how are they embedded and compared?\n\n4. In the tables of numerical results, the meaning of the reported percentages is not explained. Do they represent accuracy, relative error, or improvement over a baseline?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "The problem is challenging, and the numerical results appear promising.", "weaknesses": "There is no comparison with classical non–learning methods. For example, in the eigenvalue problem, when approximating the smallest eigenvalue, conventional numerical approaches are already highly efficient and accurate. Without such baselines, it is difficult to assess whether the proposed learning framework offers any practical advantage.\n\nThe presentation is weak, making the article difficult to read. Beyond apparent typos (see below), the referee identifies the following concerns:\n\n   1. Above Eq. (1), both matrices and parametric models are denoted by $W$, which is confusing.\n\n   2. Above Eq. (1), in the general setting, it would help to specify $V(r)$ explicitly for each example discussed later.\n\n   3. In the approximate eigenspaces example, $U$ and $E$ should be defined. On line 113, the object described is not a subspace.\n\n   4. On line 141, $V$ should be a matrix rather than a subspace based on context. If $V$ were a subspace, $S(V)$ would be ill-defined.\n\n   5. In Theorem 3, the mappings $F$ and $G$ should be indexed as $F_k$ and $G_k$, since they are $k$-dependent. The quantities $\\#_F(k, D)$ and $\\#_G(k, D)$ are not defined.\n\n   6. The sentence “As the main neural network architecture we used FFNO Tran et al. (2021), a descendant of FNO Li et al. (2020), and performed extensive grid searches for all experiments” is unclear; its relevance and meaning should be clarified.\n\n\nNumerous grammar issues:\n\n   1. Line 107: “to approximates” → “to approximate”\n   2. Line 112: “is solver repeatedly” → unclear phrasing\n   3. Line 206: “we allowed” → revise for tense consistency\n   4. Line 261: “Theorem 3 suggest” → “Theorem 3 suggests”\n   5. Line 304: missing punctuation\n\n   …and many others throughout the manuscript.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jEDi62EKmP", "reviewer": "Reviewer_MfjB", "summary": "This work proposes a reliable approach for subspace regression in ROM problems. To address the difficulty of recovering the exact number of subspace bases, the authors introduce a relaxation named subspace embedding that learns a higher-dimensional subspace instead, intended to sufficiently encompass the target lower-dimensional subspace and improve both accuracy and robustness at the expense of some gap between the dimensions of optimal and learned subspaces. This new formulation is evaluated in several applications and compared to different baselines to highlight its pros and cons.", "questions": "1. How is \"relative error\" defined? $W_\\theta$ and $V$ should have different ranks ($r$ and $k$), right? \n   * Maybe related to the Frobenius norm in Theorem 2?\n   * Is the definition of \"relative error\" the same for different models in Figure 2?\n2. Will subspace embedding perform worse as the ambient dimension $n$ increases? I assume that as $n$ increases, the low-dim linear subspaces will be more likely orthogonal with one another and make subspace embedding trickier. Is this intuition reasonable?\n   * In other words, what would Figure 1 approximately look like if $n$ is the third dimension apart from $r$ and error?\n3. How impactful is the gap $r-k$ to different ROM problems when the predicted larger subspace is used? \n   * Any postprocess for trimming the predicted subspace to reduce its rank to $k$?\n4. How trainable or reusable is a subspace regressor for a collection of problems of varying dimension $n$ or domain $\\Omega$? \n   * For instance, FEA is widely used to solve PDEs over domains of complicated geometries, and its approximate solution often takes the form on Line 107. How applicable is your method considering that the number of $\\phi_i$ often changes due to either the varying nonlinearity in the solution field or the varying $\\Omega$?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "1. Well-structured paper with clear problem definition, rigorous theoretical analysis, thorough experiments and meticulous evaluation. \n2. The subspace embedding method is novel and seems easily generalizable to different ROM problems with good robustness. \n3. Different implementations of the formulation are given in Theorem 1. \n3. Both pros and cons of this method are discussed in detail, the explanations are insightful and intuitive. \n4. Several nice illustrations helping readers appreciate the theoretical analysis in the appendix.", "weaknesses": "Overall, this is a solid paper without significant flaws. Just a few minor ones. I'm willing to increase the score if they are addressed and my questions are answered. \n\n1. Seems like a $\\psi(x)$ next to $U(x)$ is missing in (3).\n2. I think a clear definition of \"relative error\" is needed. See Question 1 below.\n3. Use \\citep and \\citet properly.\n4. Some sentences could be streamlined and a few typos can be fixed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "JX5xcOnq8t", "reviewer": "Reviewer_jxox", "summary": "The problem and a method of subspace regression (using deep learning) are proposed. The problem is to predict a subspace (i.e., an element of a Grassmann manifold) from some inputs. The authors propose loss functions, use cases of neural net-based architectures, and a way to facilitate training by regressing subspaces with larger dimensionality. The numerical examples also support the utility of the proposed framework.", "questions": "What are the most relevant studies, particularly in terms of the problem setting, i.e., regression to subspaces? As the paper does not have a distinct Related Work section, it is a little difficult to follow the exact technical context.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The motivation is clearly explained with concrete examples, the proposed framework is technically sound, and the numerical examples are convincing to show the approach's effectiveness.\n\nI assume the problem setting itself is something new, but as I am not an expert in the area, I cannot fully assess the novelty in the particular context of subspace regression.", "weaknesses": "I basically don't see major issues.\n\nAlthough the paper is carefully written, a kind of schematic drawing of (an instance of) the problem and the method will be highly appreciated to facilitate smoother understanding.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)