PMaF: Deep Declarative Layers for Principal Matrix Features

Zhiwei Xu 1 Hao Wang† 1 Yanbin Liu 1 Stephen Gould 1

Abstract
We explore two differentiable deep declarative
layers, namely least squares on sphere (LESS) and
implicit eigen decomposition (IED), for learning
principal matrix features (PMaF). It can be used
to represent data features with a low-dimensional
vector containing dominant information from ahigh-dimensional matrix. We first solve the problems with iterative optimization in the forward
pass and then backpropagate the solution for
implicit gradients under a bi-level optimization
framework. Particularly, adaptive descent steps
with the backtracking line search method anddescent decay in the tangent space are studied
to improve the forward pass efficiency of LESS.
Meanwhile, exploited data structures are used togreatly reduce the computational complexity inthe backward pass of LESS and IED. Empirically,
we demonstrate the superiority of our layers over
the off-the-shelf baselines by comparing the solution optimality and computational requirements.

Figure 1. With BLS and TWD, LESS converges faster and better.
The moving path starts from the initial value (red-dot) to the optimal (black-star). More results are in the Appendix.

for PMaF and study two deep layers, namely least squares
on sphere (LESS) and implicit eigen decomposition (IED),
that are superior to off-the-shelf SciPy (nondifferentiable)(Virtanen et al., 2020) and PyTorch (Paszke et al., 2019)baselines in the effectiveness and/or efficiency of optimization and differentiability for end-to-end learning.
In LESS, the proposed adaptive gradient descent steps on
the tangent plane greatly reduce the number of iterations,
see Fig. 1; in IED, the alternatives, power iteration (PI)(von Mises & Pollaczek-Geiringer, 1929) and simultaneous
iteration (SI) (McCormick & Noe, 1977), achieve better
solutions for non-negative symmetric and nonsymmetricmatrices than the baseline. Meanwhile, we use implicit
differentiation methods, mainly deep declarative network(DDN) (Gould et al., 2021) in this work, with exploited matrix structures (Gould et al., 2022) to dramatically reduce thecomputational requirements. Comprehensive experiments
are provided in the Appendix.

1. Introduction
Principal matrix feature (PMaF)‡ in this work refers to asingle vector summarising a data matrix. It can be usedin deep feature representation or learning that is typical invarious areas, such as image analysis (Xu et al., 2014; MelasKyriazi et al., 2022), natural language processing (Younget al., 2017), weather prediction (Malakar et al., 2021), andso on. It adapts learned features for downstream tasks andstudies matrix structures for fine-grained features with suchas a high sparsity or a low dimension (Ranzato et al., 2007;
Liu & Yan, 2011; Robles-Kelly, 2016; Liu et al., 2017).

2. Deep Declarative Layers
2.1. Least Squares on Sphere
Given A ∈ Rm×n and b ∈ Rm , the least squares problem
with solution constrained on a unit sphere is defined as

In this work, we mainly focus on two optimization problems

1
2∥Au − b∥ ,22subject to ∥u∥ = 1 ,minimizeu∈Rn f(A, b, u) ≜†Parts of this work were completed when Hao Wang pursued
a master’s degree at the Australian National University. ‡ PMaF
for the IED problem is equivalent to principal component analysis(Mackiewicz & Ratajczak, 1993) for a covariance matrix. 1 Schoolof Computing, CECC, ANU, Canberra, Australia. Correspondence
to: Zhiwei Xu <zhiwei.xu@anu.edu.au>.

(1)

where ∥·∥ is the ℓ2 -norm. For the notation simplicity, we
denote f(A, b, u) by using f(u). Since u is constrainedon the sphere, the closed-form of the optimal y = A−1 b
no longer holds. The projected gradient descent (PGD)
method is used to decrease the energy in the gradient descentdirection while guaranteeing the solution feasibility.

Published at the Differentiable Almost Everything Workshop of the
40 th International Conference on Machine Learning, Honolulu,
Hawaii, USA. July 2023. Copyright 2023 by the author(s).

1PMaF: Deep Declarative Layers for Principal Matrix Features

We first describe the vanilla PGD, PGD with direction
weight decay, and PGD projected onto the Riemannian manifold (directly on the constraint sphere). For monotonic
energy convergence, either the backtracking line searchmethod or a simple yet effective step decay in the tangent
space is used for fast optimal solution search.

iteration is crucial, requiring a suitable descent step. For
this to hold, we refer to the backtracking line search method(Boyd & Vandenberghe, 2004) and apply the first-order
Lagrangian form of Eq. (1)f(ut + η∆ut ) ≤ f(ut ) + αη∇⊤ f(ut )∆ut ,(7)1) Projected Gradient Descent (PGD). We refer to the Appendix for the details of the widely used PGD method (Akilov & Kantorovich, 1982; Lemarechal, 2012).

where constant α ∈ (0, 0.5) is for the maximum energy decrease at the (t + 1)th iteration for the descent monotonicity.
Otherwise, the descent step would be decayed by η ← βηwith β ∈ (0, 1) to avoid surpassing the optimal solution.

2) Direction Weight (DW). Since the descent step η in PGD
needs to decrease for a fine search when ut approaches the
optimal solution, reducing η with weight wt adaptive to thedescent direction is desirable. Given f(u) in Eq. (1), thecenter u0 = (A⊤ A)−1 A⊤ b indicates if the least squares
problem is an inner (non-convex) or outer (convex) equalityconstrained problem and also the direction from the constraint center to the solution.

5) Tangent Weight Decay (TWD). Alternatively, since alarge descent step tackles the monotonic energy decrease in
the update around the optimal solution, it could be unable
to converge on the manifold. Hence, the descent step isdecayed in the tangent space on the same side of the descentdirection with decay rate β, that is η ← βη, when
Sc (ProjRM (−∇f(ut )) , ProjRM (−∇f(ut+1 ))) < 0 ,(8)with
ut+1 = ut + ηProjRM (−∇f(ut )) .
(9)

With the optimal unconstrained solution, u0 is parallel to thedescent direction dt . Hence, we define a direction weight
for each descent step as
wt = 1 − Sc (dt , u0 ) ,

(2)

(

(3)

where
dt =

−∇f(ut )∇f(ut )

if ∥u0 ∥ ≥ 1otherwise

In short, Eq. (8) indicates that the solution update at the (t +1)th iteration causes a “reverse” descent direction, usuallyin (90, 180] degrees, around the optimal solution, and thus,a smaller step is preferred.
2.2. Implicit Eigen Decomposition⊤

and Sc (a, b) = a b/ (∥a∥ ∥b∥) is the cosine similarity
to measure the direction homotropy. When ut is optimal,
wt = 0 such that wt η = 0 terminates the update of ut .

Generally, the eigen decomposition problem can be formulated by solving the optimization f : A ∈ Rm×m →{λ, u} with λ as a vector of n largest eigenvalues andu ∈ Rm×n as the corresponding eigenvectors,
minimizeu∈Rm×n f(A, u) ≜ −tr u⊤ Au ,(10)subject to h (u) ≜ u⊤ u = In ,3) PGD on Riemannian Manifold (RM). As the solution
is constrained on a sphere, using Riemannian geodesic distance is promising to learn better feature distribution thanusing the Euclidean distance (Barachant et al., 2010; Wanget al., 2017). Therefore, we define the solution projection
onto the Riemannian manifold (Boumal, 2020) as
ProjRM (−∇f(ut )) = In − ut u⊤t (−∇f(ut )) , (4)where tr() is the trace function over all n eigenvalues. The
optimal solution in Eq. (11) satisfies Eq. (12) as
y = argminu∈Rm×n f(A, u) ,Ayi = λi yi , ∀i ∈ N = {1, ..., n} .

where In is an n × n identity matrix. The solution update
at (t + 1) follows
ut+1 = ProjSph (ut + ηProjRM (−∇f(ut ))) ,

(11)
(12)

The principal matrix component refers to the eigenvector
associated with the largest eigenvalue, for which we use thepower iteration algorithm and the simultaneous iteration algorithm. For a complete analysis, solution update formulas
are provided below.

(5)

where η = 1 or can be derived from ∇f(ut+1 ) = 0 as ⊤⊤(Aut − b) A In − ut u⊤t A (Aut − b)
η=
. (6)
2
⊤
A(In − ut u⊤t )A (Aut − b)

Solver 1. Power Iteration (PI). Given a randomly initialized
eigenvector u0 , the solution update at time (t + 1) follows
ut+1 = Aut / ∥Aut ∥ ,See the Appendix for the Riemannian manifold projection.

(13)

and terminates (upon the convergence or the maximum iteration) at t = K for the principal eigenvector and eigenvalue,

4) Backtracking Line Search (BLS). While the direction
weight method causes many iterations if the solution isfar from optimal, a guarantee of energy reduction in every

y = uK2and λ = y⊤ Ay .

(14)

PMaF: Deep Declarative Layers for Principal Matrix Features

Solver 2. Simultaneous Iteration (SI). QR decomposition
is required for iterative updates of the input and the solution.
The initial input x0 = A and the updates at time t follow{Qt , Rt } = QR(xt )and

xt+1 = xt Qt .

Proposition 3.1. (Exploited Hessian structure for LESS)
Rather than applying Jacobian and Hessian operations, anaccumulation by parts approach with the exploited structureof B greatly improves the implementation efficiency. Recall
that B ∈ Rn×(m×n) , A ∈ Rm×n , and y ∈ Rn . Set indicesi, j ∈ N = {1, ..., n} for the dimensions related to n suchthat B = [Bij ], A = [Ai ], and y = [yj ]. Then,(15)

The principal eigenvector is the component of QK corresponding to the largest eigenvalue λ = max(RK ).
Solution consistency for effective backpropagation. The
eigen decomposition problem defined in Eq. (10) has two
optimal solutions with reverse directions due to the quadratics. This could happen either in the intermediately updated
solution ut (either Eq. (13) or Eq. (15)) or in the optimal
solution y in each learning epoch for the same data sample.
To alleviate the ineffectiveness of these updates, one canapply a reference direction r such that the optimal solution
y updates in the same direction of r as

Proof Sketch. See the Appendix for details.

Historical : ut ← V (ut , ut−1 ) ut ,(16)
Hard-coded : y ← V (y, r) y ,(17)where V(a, b) = Sign a⊤ b if (a ̸⊥ b) and otherwise 1
and Sign() calculates the sign value of a scalar. If ut is foreigenvectors of multiple eigenvalues, only the diagonals ofu⊤t ut−1 are used for the sign values.

1
β = − y⊤ Ay + A⊤ y ∈ R .
2

Bij = Ai yj⊤ ,Bii ← Bii + (Ay − b),

(22)

(27)

Biij ← Biij − yj ,

(28)
(29)

∀i, j ∈ M .

(30)
(31)
(32)

Proof Sketch. This can be easily proved by finding thestructure of B from an example.
□
Nevertheless, B requires (4m3 /10242 )MB memory inthe single-precision floating-point format, for instance,
4,096MB for m = 1024. With the exploited Hessian structure of B and the vector-Jacobian product for ∇X L, however, B can be avoided by using merely y.

1) LESS. The components of K in Eq. (18) are

H = A⊤ A − 2βIn ∈ Rn×n ,1
β = y⊤ A⊤ (Ay − b) ∈ R .
2(26)

B = ∇2XY f(A, y) ∈ Rm×(m×m) ,
H = − A + A⊤ − 2βIm ∈ Rm×m ,Bijk = 0, ∀i, j, k ∈ M ,Biji ← Biji − yj , ∀i, j ∈ M ,

(19)

(21)

A = 2y⊤ ∈ R1×m ,Proposition 3.3. (Exploited Hessian structure for IED)
Recall that B ∈ Rm×(m×m) and y ∈ Rm . Set indicesi, j, k ∈ M = {1, ..., m} for the dimensions related to msuch that B = [Bijk ] and y = [yi ]. Then,

With the iterative optimization as the forward pass, the backward pass of y to A in Eq. (1) and Eq. (10) is required.
Backtracking the forward pass, however, is inefficient andsometimes infeasible due to the discrete solution. Hence,
we use the deep declarative network method (Gould et al.,2021) to efficiently calculate ∇X L with ∇X y, where L is
the loss from the upper problem in the bi-level optimization
and X and Y indicate the input variable (can be multiple)
and the optimization solution respectively. It follows


−1
−1
−1
−1
A − In H , (18)K = ∇Y L H A⊤ AH A⊤

B = ∇2XY f(A, b, y) ∈ Rn×(m×n) ,□Proof. Since ∇Y f(A, y) = β∇Y h(A, y) in Eq. (17) of(Gould et al., 2021) for the minima, ∇Y f(A, y) = −(Ay+A⊤ y)⊤ and ∇Y h(A, y) = 2y⊤ for the problem defined in
Eq. (10), and Ay = λy, one has −(λy + A⊤ y) = 2βy. If
A⊤ = A, then A⊤ y = Ay = λy, and thus, β = −λ.

3.1. Deep Declarative Networks based Gradients

(20)

(25)

Proposition 3.2. If the input A of the objective function
Eq. (10) is symmetric, the Lagrange multiplier equals thenegative largest eigenvalue, that is, β = −λ.

Given the optimal solutions, their gradients over the input
entries enable end-to-end learning. Without unrolling theiterations of the solvers, we use a single-step method with
implicit differentiation (Gould et al., 2021), fixed-point theorem (for IED), and exploited structures (Gould et al., 2022).

A = 2y⊤ ∈ R1×n ,∀i ∈ N .

(24)2) IED. Similarly, for IED, one has

3. Implicit Differentiation∇X L = KB .

∀i, j ∈ N ,Proposition 3.4. (Exploited gradient structure for IED) B
can be split into y dependencies for ∇X L ∈ Rm×m as∇X L = −Ky⊤ − yK⊤ .

(23)

Proof. See the Appendix for details.
3(33)

PMaF: Deep Declarative Layers for Principal Matrix Features

3.2. Implicit Function Theorem based Gradients

Table 1. With 1,000 random data from N (0, 1). “In” and “Out”
refer to the number of failed cases inner and outer of the sphererespectively. “Imp.” is the number of cases with FPD ≤ the SciPy
baseline. See the Appendix for more results.
Size 2×2
Size 64×32
MethodIn↓ Out↓ Imp.↑ MRE↓ In↓ Out↓ Imp.↑ MRE↓
PGD
172 353 389 10.84 548 440
52 0.15+RM
0 118 844 4.30 30
0 980 0.00
+RM+BLSOne 0 14 938 0.14 6
0 994 0.00
+RM+TWD00 806 0.00 29
0 954 0.00

The implicit function theorem based gradients only apply toeigen decomposition. By using PI for the largest eigenvalue(to be positive), the implicit function isf(A, ut , ut+1 ) = ut+1 − Aut / ∥Aut ∥(34)

with the input matrix A ∈ Rm×m and the solution u ∈Rm×n . Upon the convergence, y = ut+1 = ut , and thus,f(A, y) = y − Ay/ ∥Ay∥ .

(35)

Applying the implicit function theorem to f(A, y) achieves∇X f(A, y) + ∇Y f(A, y)∇X y = 0 ,∇X y = − (∇Y f(A, y))−1∇X f(A, y) .

Table 2. Backward speedup with exploited structures. Numbers
are averaged over 100 samples. “PGD+RM+TWD” is used forLESS. “Speedup” is (“AutoDiff”-“LESS”)/“LESS”.
Size 256 ×8×64
Size 256×16×128
Time (s) Memory (MB) Time (s) Memory (MB)AutoDiff
4.44
65.22
8.64
519.32
LESS
0.01
49.20
0.02
325.70Speedup ×444↑
×0.33↑
×432↑
×0.59↑(36)(37)

For the notation consistency with Sec. 3.1, we denote H =∇Y f(A, y) ∈ Rm×m , B = ∇X f(A, y) ∈ Rm×(m×m) ,and ∇X y = −H−1 B ∈ Rm×(m×m) . Note that ∇X y isin the same form of Eq. (19) as an unconstrained problem.
Then, ∇Y L ∈ Rm follows−1∇X L = ∇Y L∇X y = −∇Y LH B ,
H = Im − Im − yy⊤ A/λ ,
B = − Im − yy⊤ y⊤ /λ .

(38)
(39)
(40)

Different from using multiple eigenvalues for gradients(Magnus & Neudecker, 2007), only the largest one is used.
Again, it is strongly encouraged to use the vector-Jacobian
product to calculate ∇X L over B for memory reduction.

Figure 2. Precision evaluation with symmetric A in float32. See
the Appendix for nonsymmetric A and simulation on ResNet50.
Table 3. Efficiency evaluation with symmetric A and solution size
5 × 512. “J”: autodiff Jacobian without exploited structures; “E”:
our differentiation with exploited structures. See the Appendix forresults with nonsymmetric A and more sizes.
GPU Time (s)GPU Memory (MB)MethodForward Backward Forward Backward
AutoDiff
0.07580.0005
20.56
24.99
PI-DDN-J
0.00690.0436
5.55
2570.06
0.00700.69425.55
5125.00
PI-IFT-J
PI-DDN-E
0.00690.0114
5.55
16.01
PI-IFT-E
0.00690.0043
5.55
31.01
SI-DDN-J
1.06230.041932.02
2570.06
1.06240.712032.02
5125.00
SI-IFT-J
SI-DDN-E
1.06220.0113
32.02
16.01
SI-IFT-E
1.06250.004232.02
31.01

4. Experiments
Experiments are on a 3090 GPU and PyTorch 1.12.0. The
code is at https://github.com/anucvml/ddn.git. See the Appendix for evaluation standards including fixed-point distance (FPD) in Eq. (41), eigen distance in Eq. (43), andmean relative error (MRE) in Eq. (44). We highlight ours .
4.1. Evaluation on LESS. We show the effectiveness in
Fig. 1 and Table 1 and efficiency in Table 2. Fig. 1 showsa case where ours can achieve fewer iterative steps when
optimizing Eq. (1). In Table 1, ours with BLSOne and
TWD achieve comparable low MRE with the SciPy solver(Virtanen et al., 2020) (supports only CPU and not differentiable) and outperform the vanilla PGD and PGD+RM.
With exploited structures including Prop. 3.1, ours is fasterand more memory efficient than without those structures(“AutoDiff”). More results are in the Appendix.

5. Conclusion
We explore two deep layers, LESS and IED, for PMaF thatcan be used in end-to-end learning. Either optimization improvements or alternatives enhance the optimization effectiveness and/or efficiency. We also enable the differentiation
feasibility and largely reduce computational requirements
using DDN and IFT based gradients with exploited structures. Extensive experiments show the benefits of our layers
over the baseline from SciPy optimizer or PyTorch function.

4.2. Evaluation on IED. We show the effectiveness in Fig. 2and efficiency in Table 3. In Fig. 2, PI and SI achieve lower
eigen distance and FPD than PyTorch eigh() (“AutoDiff”)(Paszke et al., 2019). With exploited structures in Table 3,ours achieve higher speed or less memory without ignoring
their effectiveness. More results are in the Appendix.
4

PMaF: Deep Declarative Layers for Principal Matrix Features

References

Malakar, S., Goswami, S., Ganguli, B., Chakrabarti, A.,
Roy, S. S., Boopathi, K., and Rangaraj, A. G. A novel
feature representation for prediction of global horizontal irradiance using a bidirectional model. In Machine
Learning and Knowledge Extraction, 2021.

Akilov, G. P. and Kantorovich, L. V. Functional analysis(2nd edition). In Pergamon Press, 1982.
Barachant, A., Bonnet, S., Congedo, M., and Jutten, C.
Riemannian geometry applied to BCI classification. In
nternational Conference on Latent Variable Analysisand
Signal Separation, 2010.

McCormick, S. F. and Noe, T. D. Simultaneous iteration forthe matrix eigenvalue problem. In Linear Algebra and its
Applications, 1977.

Boumal, N. An introduction to optimization on smooth manifolds. http://www.nicolasboumal.net/book, Aug 2020.

Melas-Kyriazi, L., Rupprecht, C., Laina, I., and Vedaldi, A.
Deep spectral methods: A surprisingly strong baseline
for unsupervised semantic segmentation and localization.
In IEEE International Conference on Computer Vision,
2022.

Boyd, S. and Vandenberghe, L. Convex optimization. Cambridge University Press, 2004.

Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J.,
Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga,
L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang,
L., Bai, J., and Chintala, S. PyTorch: An imperative style,high-performance deep learning library. In Advancesin Neural Information Processing Systems, 2019.
URL https://pytorch.org/docs/stable/
generated/torch.linalg.eigh.html.

Gould, S., Hartley, R., and Campbell, D. Deep declarative
networks. In IEEE Transactions on Pattern Analysis andMachine Intelligence, 2021.
Gould, S., Campbell, D., Ben-Shabat, I., Koneputugodage,
C. H., and Xu, Z. Exploiting problem structure in deep
declarative networks: Two case studies. In AAAI Workshop on Optimal Transport and Structured Data Modeling, 2022.

Ranzato, M., Boureau, Y. L., and LeCun, Y. Sparse featurelearning for deep belief networks. In Conference on
Neural Information Processing Systems, 2007.

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M.,
Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane,A., del Rı́o, J. F., Wiebe, M., Peterson, P., GérardMarchant, P., Sheppard, K., Reddy, T., Weckesser, W.,
Abbasi, H., Gohlke, C., and Oliphant, T. E. Array programming with NumPy. In Springer Science and Business Media LLC, volume 585, pp. 357–362, September
2020. doi: 10.1038/s41586-020-2649-2. URL https://doi.org/10.1038/s41586-020-2649-2.

Robles-Kelly, A. Least-squares regression with unitary constraints for network behaviour classifications. In Joint
IAPR International Workshops on Statistical Techniques
in Pattern Recognition and Structural and Syntactic Pattern Recognition, 2016.
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland,M., Reddy, T., Cournapeau, D., Burovski, E., Peterson,
P., Weckesser, W., Bright, J., van der Walt, S. J.,
Brett, M., Wilson, J., Millman, K. J., Mayorov, N.,
Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey,
C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas,J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen,
I., Quintero, E. A., Harris, C. R., Archibald, A. M.,
Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and
SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. In Nature
Methods, 2020. URL https://docs.scipy.org/
doc/scipy/reference/generated/scipy.
optimize.minimize.html.

Lemarechal, C. Cauchy and the gradient method. In Doc
Math Extra, 2012.
Liu, G. and Yan, S. Latent low-rank representation for
subspace segmentationand feature extraction. In IEEE
International Conference on Computer Vision, 2011.
Liu, J., Wu, Z., Xiao, Z., and Yang, J. Classification of
hyperspectral images using kernel fully constrained leastsquares. In International Journal of Geo-Information,
2017.
Mackiewicz, A. and Ratajczak, W. Principal components
analysis (PCA). In Computers and Geosciences, 1993.

von Mises, R. Mathematical theory of probability and statistics. In New York: Academic Press, 1964.

Magnus, J. R. and Neudecker, H. Matrix differential calculus with applications in statistics and econometrics (3rdedition). In John Wiley and Sons Ltd, 2007.

von Mises, R. and Pollaczek-Geiringer, H. Praktische verfahren der gleichungsauflosung. In Zeitschrift fur Angewandte Mathematik und Mechanik, 1929.
5

PMaF: Deep Declarative Layers for Principal Matrix Features

Wang, Y., Li, W., Dai, D., and Gool, L. V. Deep domain
adaptation by geodesic distance minimization. In IEEE
International Conference on Computer Vision, 2017.
Xu, Y., Mo, T., Feng, Q., Zhong, P., Lai, M., and Chang,
E. I.-C. Deep learning of feature representation withmultiple instancelearning for medical image analysis. In
IEEE International Conference on Acoustics, Speech andSignal Processing, 2014.
Young, T., Hazarika, D., Poria, S., and Cambria, E. Recent
trends in deep learning based natural language processing.
In arXiv:1708.02709, 2017.

6

PMaF: Deep Declarative Layers for Principal Matrix Features

Appendix
A. Evaluation Standards
1) Fixed Point Distance (FPD). This measures the distance of the estimated solution to the minimum of the objectivefunction, where the minimum objective value is unnecessary to be zero given problem constraints.
LESS :

F(y) = ∥Ay − b∥ ,IED :F(y) = y −

Ay
∥Ay∥

(41)
.

(42)

2) Eigen Distance. This is intuitive from Eq. (12) to evaluate the eigen solution,
D(y) = ∥Ay − λy∥ .

(43)3) Mean Relative Error (MRE). Given a reference (generally ground truth) yr and k th sample, the mean relative error isK

1 X F(yk ) − F(ykr )× 100% ,MRE(y, y ) =KF(ykr )
r

(44)

k=1

which can be negative when the estimation y outperforms the reference yr .
B. Projected Gradient Descent and Riemannian Manifold Projection

Figure 3. Descent direction projected onto a Riemannian manifold for the sphere constraints. −∇f(xt ) is the PGD direction, G =ProjRM (−∇f(xt )) is the descent direction on the Riemannian manifold, and xt is on the Riemannian manifold.

B.1. P ROJECTED G RADIENT D ESCENT
For ease of optimization, the solution is initialized with the least squares solution without constraints as
−1 ⊤
A b,
u 0 = A⊤ A
u0 ← ProjSph (u0 ) .

(45)
(46)

Then, the solution is iteratively updated in the gradient descent direction ∆ut at tth iteration by following
ut+1 = ut + η∆ut ,(47)

ut+1 ← ProjSph (ut+1 )

(48)

with
∆ut = −∇f(ut ) ,η=A⊤ (Aut − b)2
2 ,∥AA⊤ (Aut − b)∥(49)

where η is derived by setting ∇f(η) ≡ ∇f(ut+1 ) = 0 and ProjSph (u) = u/ ∥u∥ is the projection onto the unit sphere.
The solution update terminates when the objective value decrease is no greater than a predefined tolerance, 1.0e−7 in ourcase, or when it reaches the maximum number of iterations, 100 in our case.
7

PMaF: Deep Declarative Layers for Principal Matrix Features

B.2. R IEMANNIAN M ANIFOLD P ROJECTION
Since the solution is constrained on the sphere, a simple projection onto the sphere via gradient descent is inefficientconsidering the inactive orthogonal descent component. In Fig. 3, considering only the descent component on the tangent
plane benefits the solution optimality and optimization process, comparing “PGD” and “PGD+RM” in Table 4.
−xt∇f(xt ) ∥x
−xt
−xt
t∥
∥∇f(xt )∥ cosθ =∥∇f(xt )∥
P=
∥xt ∥
∥xt ∥
∥∇f(x )∥ xt
t

then

=

∥xt ∥

xt x⊤
t

2 ∇f(xt ) ,

∥xt ∥

xt x⊤
t

G = −∇f(xt ) + P = −∇f(xt ) +(50)

2 ∇f(xt ) .

(51)(−∇f(xt )) ,

(52)

∥xt ∥

Hence,
ProjRM (−∇f(xt )) =

In −

xt x⊤
t
2

∥xt ∥

!

where xt ∈ Rn , ∇f(xt ) ∈ Rn , In is an n × n identity matrix.
C. Gradient Transformation with Flipped Eigenvalues
For different orderings of multiple eigenvalues, for instance, descending-order (concerning the order of eigenvalues) versus
ascending-order eigenvalues, the gradients of their associated eigenvectors need to be flipped according to the flipping of
the eigenvalues. Recall the eigen decomposition problem defined in Eq. (10). Reserve all of the eigenvalues of A, that isn = m, a flipping function is defined as g(x) : Rm×m → Rm×m , that is 
0 ··· 10 ··· 1

 

(53)
g(x) =  ... . . . ...  x  ... . . . ...  = I∗m xI∗m ,
1

···

1

0

···

0

where (I∗n x) flips x in the vertical direction and (xI∗n ) in the horizontal direction. For the notation simplicity, we denote
g(x) = xflip . Then,
flip flip
Au = Λu = Iflip
m Λ Im u ,

Iflip
mA

flip flip
Iflip
m Au = Λ Im u ,flip flip
flip flip flip
Iflip
m u Im = Λ u I m ,flip flip
flip flip

A u

=Λ u

(54)
(55)
(56)

,

(57)

where the input A ∈ Rm×m , eigenvectors u ∈ Rm×m , and eigenvalue matrix Λ ∈ Rm×m . It indicates that to obtain
the reverse-order eigenvalues, A needs to be flipped, which also leads to the flipping of u. Therefore, to compare thedescending-order eigen decomposition with an ascending-order one, both A and u need to be flipped.
flip
For the gradient of loss L over the flipped eigenvectors uflip , that is ∇uflip L, since uflip = Iflip
n uIn ,
flip
∇uflip L = Iflipn ∇uflip LIn = (∇u L)

flip

.

(58)

In sum, for the reverse-order eigenvalues, gradients of their associated eigenvectors need to be flipped in both the horizontal
and vertical directions, that is to apply Eq. (53), when the eigenvalue order (or equivalently the eigenvector order) is flipped.
D. Proof and Demo Code for Proposition 3.1
We explore the matrix structure by taking an example with m = 2 and n = 3, that is A, b, and y as
 

 

y0
b0
a00 a01 a02, y =  y1  .
, b=
A=
b1
a10 a11 a12
y2
8

(59)

PMaF: Deep Declarative Layers for Principal Matrix Features

The Jacobian vector and Hessian matrices of the objective function in Eq. (1) over the solution y and the Hessian matrix
over entries A and b are∇Y f(A, b, y) = A⊤ (Ay − b)
 2

 
a00 b0 + a10 b1(a00 + a210 )y0 + (a00 a01 + a10 a11 )y1 + (a00 a02 + a10 a12 )y2= (a01 a00 + a11 a10 )y0 + (a201 + a211 )y1 + (a01 a02 + a11 a12 )y2  − a01 b0 + a11 b1  ,a02 b0 + a12 b1(a02 a00 + a12 a10 )y0 + (a02 a01 + a12 a11 y1 + (a202 + a212 )y2 )
a00 a10∇2bY f(A, b, y) = − a01 a11  = −A⊤ ,a02 a12  
b0 0 02a00 y0 + a01 y1 + a02 y2 a00 y1 a00 y2
 2a10 y0 + a11 y1 + a12 y2 a10 y1 a10 y2   b1 0 0 


 


 

 





0
b
0
a
y
a
y
+
2a
y
+
a
x
a
y
0
01
0
00
0
01
1
02
2
01
2
 

∇2AY f(A, b, y) =  a11 y0 a10 y0 + 2a11 y1 + a12 x2 a11 y2  −  0 b1 0 

 

 



  
 a02 y0 a02 y1 a00 y0 + a01 y1 + 2a02 y2   0 0 b0 
0 0 b1
a12 y0 a12 y1 a10 y0 + a11 y1 + 2a12 y2
 

 

a00 y0 a00 y1 a00 y2
a00 y0 + a01 y1 + a02 y2 0 0
b0 0 0 a10 y0 a10 y1 a10 y2   a10 y0 + a11 y1 + a12 y2 0 0   b1 0 0 

 

 


 
 

 
 

 a01 y0 a01 y1 a01 y2   0 a00 y0 + a01 y1 + a02 y2 0   0 b0 0 
 

 
=
 a11 y0 a11 y1 a11 y2  +  0 a10 y0 + a11 y1 + a12 y2 0  −  0 b1 0 
 


 
 

 


 

 
 a02 y0 a02 y1 a02 y2   0 0 a00 y0 + a01 y1 + a02 y2   0 0 b0 
a12 y0 a12 y1 a12 y2
0 0 a10 y0 + a11 y1 + a12 y2
0 0 b1  
a00 ⊤ a10 y 


 

 

  
Ay 00
b 0 0
 a01 ⊤ 
 
 

=
 a11 y  + 0 Ay 0 − 0 b 0


0
0
Ay
0
0
b


  
 a02 ⊤ 
y
a12

 

A0 y ⊤
Ay − b00
Ay − b0  ,=  A1 y ⊤  +  0⊤00
Ay − b
A2 y

(60)

(61)

(62)

where the left term can be written as Ai for all i ∈ N = {1, ..., n} and the right term assigns (Ay − b) to the diagonal of∇2AY f(A, b, y) which is B in Eq. (21). One can represent it as
Bij = Ai yj⊤ ,∀i, j ∈ N ,Bii ← Bii + (Ay − b),∀i ∈ N .

(63)
(64)

One can also obtain the same structure from different examples. It is clear that forming B using A, b, and y in this way
avoids calculating the zero vectors and matrices, and thus, becoming more efficient than using Jacobian to calculate every
element of such a spare matrix. To this end, we arrive at the result in Prop. 3.1. □9PMaF: Deep Declarative Layers for Principal Matrix Features

Meanwhile, we show an implementation example comparing the automatic Jacobian (without exploited structures) and ourswith the exploited structure in Listing 1.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

Listing 1.
An implementation example of exploited Hessian structure for LESS.
import torchfrom torch.autograd import gradfrom torch.autograd.functional import jacobian as J
def obj_fnc(A, u, b):
Au_b = torch.einsum(’bmn,bn->bm’, A, u) - b
loss = 0.5 * torch.einsum(’bm,bm->b’, Au_b, Au_b)
return loss.sum()# ==== Autogradientdef dfdy_auto_fnc(A, u, b):
return grad(obj_fnc(A, u, b), (u), create_graph=True)[0]def dffdAy_auto_fnc(A, u, b):m, n = A.shape[1:3]
D = []
with torch.enable_grad():for A_p, u_p, b_p in zip(A, u, b):
A_p = A_p.view(1, m, n)u_p = u_p.view(1, n)b_p = b_p.view(1, m)
D.append(J(lambda A: dfdy_auto_fnc(A, u_p, b_p), (A_p)))
D = torch.cat(D, dim=0)
return D# ==== Structure exploited
def dfdy_fnc(A, u, b):
Au_b = torch.einsum(’bmn,bn->bm’, A, u) - b
return torch.einsum(’bmn,bm->bn’, A, Au_b)
def dffdAy_fnc(A, u, b):batch, m, n = A.shape
D = A.new_zeros(batch, n, m, n)with torch.no_grad():for i in range(n):D_1 = torch.einsum(’bm,bn->bmn’, A[:, :, i], u)
D_2 = torch.einsum(’bmn,bn->bm’, A, u) - b
D[:, i] = D_1
D[:, i, :, i] += D_2
return D

E. Proof for Proposition 3.4
For the principal matrix features, the implicit differentiation is on the eigenvector associated with the largest eigenvalue as


−1
−1
−1
−1
−1
AH − H∇X L(y) = ∇Y L(y)∇X y = ∇Y L(y) H AT AH A⊤
B = KB ∈ Rm×m ,(65)

where K ∈ Rm is a vector given the upper-level loss L ∈ R from the bi-level problem. Then,∇X L(y) = −Ky⊤ − yK⊤ .
10

(66)

PMaF: Deep Declarative Layers for Principal Matrix Features

We find the matrix structure by taking the example of m = 2, where K = [k0 , k1 ]⊤ and y = [y0 , y1 ]⊤ . Then, 
−y0 −y1
0
0
−y0
0
−y1
0
B = B0 + B1 =+00−y0 −y1
0
−y0
0
−y1

(67)

and


−k0 y0
KB = KB0 + KB1 =−k1 y0 −k0 y1−k0 y0+−k1 y1−k0 y1 −k1 y0k = − 0 y0
−k1 y1k1 y 
y1 − 0 k 0y1
k1 = −Ky⊤ − yK⊤ ,(68)
where KB0 is the outer-product of K and (−y), both KB1 and KB0 will be reshaped to m × m matrices. One can also apply
the same rule to different examples, which will have the same structure in Eq. (66). Clearly, using y ∈ Rm for ∇X L(y) ismuch more efficient and requires much less memory than using B ∈ Rm×(m×m) . To this end, we prove Prop. 3.4. □F. More Results of LESS Optimization
In addition to Fig. 1, Fig. 4 contains the full list of the compared methods. The inputs for Fig. 1 are
0.569525 −1.254572−1.583332
A=and b =.
0.414020 0.124439−0.286124(69)

In this case, the vanilla PGD requires more iterations than our “PGD+RM+BLSOne” and “PGD+RM+TWD”. Statistically,
from Table 4, when comparing “PGD” with “PGD+RM”, all the metrics are better than “PGD”. With “+BLSOne” and“+TWD”, these can be further improved.
Table 4. Effectiveness evaluation on 1,000 random data sampled from N (0, 1) normal distribution. Bold: the best, underline: the second
best. “PGD”: projected gradient descent (maximum 100 iterations), “RM”: Riemannian manifold, “BLS”: backtracking line search(α=0.5 and β=0.8), “TWD”: tangent weight decay on Riemannian manifold (β=0.9), “DW”: direction weight, “In” and “Out”: thenumber of failed cases inner and outer of the sphere respectively, where the case is regarded as failed when the solution update reaches
the maximum 100 iterations, “Imp.”: the number of cases with energy no greater than SciPy. Problem sizes, m×n for A, are 2×2 (282inner and 718 outer), 64×32 (555 inner and 445 outer), and 1024×256 (all inner).

SciPy PGD RM
✓
✓
✓
✓
✓
✓
✓
✓
✓
✓

✓
✓
✓
✓
✓
✓
✓

BLS
2×2
64×32
1024×256
TWD DW
η=1 η=∇In↓ Out↓ Imp.↑ MRE↓ In↓ Out↓ Imp.↑ MRE↓ In↓ Out↓ Imp.↑ MRE↓
172 353 389 10.84 548 440
520.15 1000 0
1.35✓ 179 43837
4.89 474 363
63
0.07 1000 00.47
0 118 844
4.30 30
0 980
0.00
4 1000 -0.03✓
014 938
0.14
6
0 994
0.00
4 1000 -0.03
40
3.51 285 142 286
0.01 984 999 -0.03✓✓ 175 421✓
0 108 853
4.21
7
0 994
0.00
2 1000 -0.03
17
1.19 472 355
92
0.04 1000 00.42✓✓ 204 536✓00 806
0.00 290 954
0.00
4 1000 -0.03
15
3.12 466 352
900.04 1000 00.42✓✓ 204 533

11

PMaF: Deep Declarative Layers for Principal Matrix Features

Figure 4. Illustration of the iterative optimization in LESS. Each row is a sample with 6 methods. “blue dash”: the least squares function;“red solid”: the sphere constraints; “green solid”: solution updates before and after ProjRM (); “green dot-dash”: the least squares
function with the final solution; “red dot”: the initial solution; “black star”: the final solution, not always the optimal. The leastsquares function with the exact solution (green dot-dash) should be tangent to the sphere (red solid). Our “PGD+RM+BLSOne” and“PGD+RM+TWD” require much fewer iterations than the others for the optimal solution with comparable fixed point distance (FPD).

12

PMaF: Deep Declarative Layers for Principal Matrix Features

G. Evaluation of IED Precision and Computational Requirements
G.1. S YMMETRIC AND N ONSYMMETRIC DATA S AMPLING
Sampling Distributions. Our evaluated data A are sampled from the standard Gaussian distribution N (0, 1), uniform distribution in [0, 1), von Mises distribution in V(0, 1) (von Mises, 1964), and random choices from [0.0, 10.0), with absolute values, where A ← A + A⊤ is applied for symmetry. We use numpy.random.randn(), numpy.random.uniform(),
numpy.random.vonmises(0, 1), and numpy.random.choice(10.0), respectively, from the NumPy library (Harris
et al., 2020). They achieve similar numerical ranges as that from the Gaussian distribution in Figs. 5-6. Hyperparameters
required by the sampling methods are not limited to the settings in this work, such as the mean and variance in the von Mises
distribution and the value upper bound in the random choices.

(a) ED, float32(b) ED, float64(c) FPD, float32(d) FPD, float64
Figure 5. Symmetric A sampled from Gaussian distribution with activation, precision evaluation with eigen distance and FPD.

(a) ED, float32(b) ED, float64(c) FPD, float32(d) FPD, float64
Figure 6. Nonsymmetric A sampled from Gaussian distribution with activation, precision evaluation with eigen distance and FPD.

(a) Symmetric A, ED(b) Symmetric A, FPD(c) Nonsymmetric A, ED(d) Nonsymmetric A, FPD
Figure 7. Evaluation on 3 data sizes within the memory capacity as it requires a (2048 × 82 ) × m2 dimension fully-connected layer in
ResNet50, where m = {32, 64, 128}. All data is in float32. Both metrics are the less the better.

The feasibility lies in 1) the simulation of normalized outputs (such as logits or softmax probability) of neural networks with
an absolute operation, 2) the requirement of positive largest eigenvalue for IFT in Eq. (34). For nonsymmetric sampling, as
we test 10,000 random samples (each with batch size 5) under this setting, the power iteration and simultaneous iteration
algorithms can always achieve almost 0 eigen distance and 0 fixed-point distance. In contrast, PyTorch eigh()* , denoted as“AutoDiff”, has a much larger eigen distance, particularly in float32.
Simulation on ResNet50. We further generate random samples by using a modified ResNet50 with pretrained weights,
removing the adaptive average pooling layer and adjusting the last fully-connected layer to obtain m2 out channels (instead* PyTorch eigh() was designed for symmetric A and PyTorch eig() for complex and general data for nonsymmetric A, but we are 1)
not aim at complex data and 2) to avoid using eigh() for symmetric matrix but rather PI or SI by investigating these results.

13

PMaF: Deep Declarative Layers for Principal Matrix Featuresof 1,000) which will then be reshaped to a m × m matrix, followed by the absolute operation and symmetry for A if required.
It follows: random data with size (5 × 3 × 256 × 256) that is (batch × channel × height × width) sampled from N (0, 1) →ResNet() → reshape to m × m → absolute operation → (A ← A + A⊤ ) if symmetric A is required. In Figs. 7, “PI” and“SI” still achieve better eigen solutions than “AutoDiff” on both symmetric and nonsymmetric A.
G.2. A DDITIONAL E XPERIMENTS ON C OMPUTATIONAL R EQUIREMENTS
We evaluate the solution precision with eigen distance and fixed-point distance on symmetric and nonsymmetric A in Fig. 5
and Fig. 6 respectively. Additionally, the computational requirements include running time and GPU memory in Tables 5-10
for both symmetric and nonsymmetric A. Although “AutoDiff” achieves faster speed in most cases, its precision is inferior
to “PI” and “SI”, particularly in float32. For the evaluation completion, we provide all of these results.
Meanwhile, since a solver achieves the same solutions given the same inputs, it consumes the same (at least quite similar)
computational resources in the backward pass. We modulate each solver with different backward propagation methods into
individual layers and provide all the forward and backward results.
In all these tables, implicit differentiation without exploited structure (“J”) and with our exploited structure (“E”) clearly
distinguishes the benefits of exploiting the Jacobian and Hessian matrices in both the running time and memory requirements.
For those with out-of-memory issues or extensive running time, we mark the results with “-” in the tables.
Table 5. Implicit eigen decomposition (IED) layer, CPU time (s) for symmetric A. “AutoDiff”: PyTorch eigh() function; “DDN”:
deep declarative network; “IFT”: implicit function theorem; “unroll”: unrolling the forward iteration for gradients via PyTorch
autodiff mechanism; “J”: autodiff Jacobian without exploited structure; “E”: our implicit differentiation with exploited structure. Our
best suggestions are highlighted considering the overall solution precision and computational requirements.
MethodAutoDiff-unroll
AutoDiff-DDN-E
PI-unroll
PI-DDN-J
PI-IFT-J
PI-DDN-E
PI-IFT-E
SI-unroll
SI-DDN-J
SI-IFT-J
SI-DDN-E
SI-IFT-E

Forward0.00040.00040.00470.00380.00370.00370.00380.01380.01270.01270.01270.01275×32
Backward
0.00010.04910.01210.00140.02300.00050.00030.01330.00140.02300.00050.0003Forward0.00070.00080.00490.00390.00390.00390.00400.04010.03870.03830.03840.0393

5×64
Backward
0.0002
0.09420.01750.00260.04520.00070.00040.02060.00260.04520.00070.0004

Forward0.00220.00210.00510.00420.00410.00410.00410.11980.11590.1096
0.1131
0.11395×128
Backward
0.00030.19090.03440.00960.11030.00120.00070.04370.0096
0.10980.00120.0007

Forward0.00750.00800.00580.00490.00520.00480.0056
0.39270.35520.37910.3666
0.36675×256
Backward
0.00112.0692
0.09840.0462
1.02830.00250.00180.13220.04641.04310.00280.0015

Table 6. Implicit eigen decomposition (IED) layer, GPU time (s) for symmetric A. “AutoDiff”: PyTorch eigh() function; “DDN”:
deep declarative network; “IFT”: implicit function theorem; “unroll”: unrolling the forward iteration for gradients via PyTorch
autodiff mechanism; “J”: autodiff Jacobian without exploited structure; “E”: our implicit differentiation with exploited structure. Our
best suggestions are highlighted considering the overall solution precision and computational requirements.
5×32
5×64
5×128
5×256
5×512
5×1024
Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward
AutoDiff-unroll
0.00040.0003 0.0041
0.0003 0.0097
0.0003 0.0261
0.0003 0.0758
0.0005 0.2193
0.0008AutoDiff-DDN-E 0.0004
0.0784 0.0042
0.1523 0.0102
0.2946 0.0274
0.5861 0.0765
1.1695 0.2258
2.3861
PI-unroll
0.00770.0168 0.0077
0.0168 0.0079
0.0171 0.0079
0.0174 0.0082
0.0172 0.0086
0.0174
PI-DDN-J
0.00650.0028 0.0066
0.0050 0.0068
0.0096 0.0067
0.0182 0.0069
0.0436 0.0068
0.5281
PI-IFT-J
0.00650.0427 0.0065
0.0840 0.0067
0.1685 0.0068
0.3313 0.0070
0.6942
PI-DDN-E
0.00650.0065 0.0065
0.0067 0.0067
0.0033 0.0068
0.0059 0.0069
0.0114 0.0070
0.0414
PI-IFT-E
0.00650.0007 0.0066
0.0009 0.0068
0.0013 0.0068
0.0023 0.0069
0.0043 0.0070
0.0152SI-unroll
0.05910.0173 0.1831
0.0169 0.2679
0.0227 0.4631
0.0580 1.0811
0.1443 2.6506
0.4876
SI-DDN-J
0.05900.0029 0.1847
0.0051 0.2694
0.0089 0.4642
0.0173 1.0623
0.0419 2.5715
0.1566
SI-IFT-J
0.05880.0426 0.1833
0.0835 0.2677
0.1663 0.4642
0.3380 1.0624
0.7120
SI-DDN-E
0.05840.0045 0.1832
0.0021 0.2676
0.0032 0.4646
0.0058 1.0622
0.0113 2.5676
0.0396
SI-IFT-E
0.05830.0007 0.1833
0.0009 0.2676
0.0013 0.4641
0.0022 1.0625
0.0042 2.6315
0.0144
Method14

PMaF: Deep Declarative Layers for Principal Matrix Features

Table 7. Implicit eigen decomposition (IED) layer, GPU memory (MB) for symmetric A. “AutoDiff”: PyTorch eigh() function; “DDN”:
deep declarative network; “IFT”: implicit function theorem; “unroll”: unrolling the forward iteration for gradients via PyTorch
autodiff mechanism; “J”: autodiff Jacobian without exploited structure; “E”: our implicit differentiation with exploited structure. Our
best suggestions are highlighted considering the overall solution precision and computational requirements.
5×32
5×64
5×128
5×256
5×512
5×1024
Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward
AutoDiff-unroll
0.5903 0.0972 0.8247 0.3896 1.7671 1.5605 5.5298 6.2456 20.5552 24.9907 80.6060
99.9810
AutoDiff-DDN-E 0.5903 0.0884 0.8232 0.3467 1.7646 1.3789 5.5249 5.5063 20.5454 22.0112 80.5864
88.0210
PI-unroll
0.8228 0.0537 0.9810 0.2280 1.4146 0.9282 2.8506 3.7334
7.5977 14.9688 24.5918
59.9395
PI-DDN-J
0.5293 0.6743 0.5898 5.1694 0.8281 40.6440 1.7759 322.5327
5.5464 2570.0620 20.5874 20520.1206
PI-IFT-J
0.5293 1.2705 0.5898 10.0791 0.8281 80.3135 1.7759 641.2510
5.5464 5125.0010
PI-DDN-E
0.5293 0.0674 0.5898 0.2534 0.8281 1.0044 1.7759 4.0068
5.5464 16.0117 20.5874
64.0215
PI-IFT-E
0.5293 0.1230 0.5898 0.4868 0.8281 1.9409 1.7759 7.7559
5.5464 31.0107 20.5874 124.0205
SI-unroll
4.6211 0.0962 16.5869 0.3882 65.0103 1.5576 256.5151 6.2402 1022.0249 24.9805 4083.0444
99.9609
SI-DDN-J
0.7539 0.6743 1.1182 5.1694 3.1353 40.6440 9.0151 322.5327 32.0249 2570.0620 123.0444 20520.1206
SI-IFT-J
0.7539 1.2705 1.1182 10.0791 3.1353 80.3135 9.0151 641.2510 32.0249 5125.0010
SI-DDN-E
0.7539 0.0674 1.1182 0.2534 3.1353 1.0044 9.0151 4.0068 32.0249 16.0117 123.0444
64.0215
SI-IFT-E
0.7539 0.1230 1.1182 0.4868 3.1353 1.9409 9.0151 7.7559 32.0249 31.0107 123.0444 124.0205
MethodTable 8. Implicit eigen decomposition (IED) layer, CPU time (s) for nonsymmetric A. “AutoDiff”: PyTorch eigh() function; “DDN”:
deep declarative network; “IFT”: implicit function theorem; “unroll”: unrolling the forward iteration for gradients via PyTorch
autodiff mechanism; “J”: autodiff Jacobian without exploited structure; “E”: our implicit differentiation with exploited structure. Our
best suggestions are highlighted considering the overall solution precision and computational requirements.
MethodAutoDiff-unroll
AutoDiff-DDN-E
PI-unroll
PI-DDN-J
PI-IFT-J
PI-DDN-E
PI-IFT-E
SI-unroll
SI-DDN-J
SI-IFT-J
SI-DDN-E
SI-IFT-E

Forward0.00040.00040.00450.00360.00370.00360.00360.01360.01270.01250.01260.01275×32
Backward
0.0001
0.04830.01130.00140.02210.00040.00030.01310.00140.02200.00040.0003Forward0.00080.00080.00460.00370.00380.00370.00380.0396
0.03870.03800.03880.03855×64
Backward
0.0002
0.09250.01500.0026
0.04320.00060.00040.01980.00270.04280.00060.0004

Forward0.00220.00210.00480.00390.00400.00390.00390.11590.11820.11040.11620.11845×128
Backward
0.00030.19020.03060.00960.10400.00120.00060.04260.00970.10450.00120.0007

Forward0.00750.00750.00540.00450.00500.00440.0046
0.36610.36510.39770.36170.37395×256
Backward
0.0011
2.09290.09230.04630.95390.00240.00150.12980.04681.00480.00280.0016

Table 9. Implicit eigen decomposition (IED) layer, GPU time (s) for nonsymmetric A. “AutoDiff”: PyTorch eigh() function; “DDN”:
deep declarative network; “IFT”: implicit function theorem; “unroll”: unrolling the forward iteration for gradients via PyTorch
autodiff mechanism; “J”: autodiff Jacobian without exploited structure; “E”: our implicit differentiation with exploited structure. Our
best suggestions are highlighted considering the overall solution precision and computational requirements.
5×32
5×64
5×128
5×256
5×512
5×1024
Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward
AutoDiff-unroll
0.00040.0003 0.0041
0.0003 0.0097
0.0003 0.0262
0.0003 0.0759
0.0005 0.2188
0.0008AutoDiff-DDN-E 0.0004
0.0816 0.0042
0.1586 0.0102
0.3062 0.0275
0.5802 0.0766
1.1564 0.2271
2.3584
PI-unroll
0.00720.0166 0.0072
0.0165 0.0074
0.0166 0.0074
0.0166 0.0077
0.0168 0.0080
0.0170
PI-DDN-J
0.00600.0028 0.0061
0.0050 0.0063
0.0090 0.0063
0.0174 0.0064
0.0431 0.0063
0.5642
PI-IFT-J
0.0061
0.0411 0.0061
0.0811 0.0065
0.1613 0.0064
0.3218 0.0066
0.6743
PI-DDN-E
0.00610.0064 0.0061
0.0068 0.0069
0.0033 0.0063
0.0059 0.0065
0.0113 0.0067
0.0410
PI-IFT-E
0.00600.0006 0.0060
0.0008 0.0064
0.0012 0.0064
0.0021 0.0065
0.0041 0.0066
0.0151SI-unroll
0.05830.0166 0.1830
0.0163 0.2674
0.0226 0.4646
0.0582 1.0884
0.1454 2.6464
0.4886
SI-DDN-J
0.05880.0029 0.1847
0.0050 0.2692
0.0089 0.4624
0.0173 1.0584
0.0420 2.5674
0.1567
SI-IFT-J
0.05900.0403 0.1832
0.0793 0.2671
0.1581 0.4623
0.3211 1.0582
0.6744SI-DDN-E
0.05840.0045 0.1833
0.0021 0.2671
0.0031 0.4625
0.0057 1.0583
0.0112 2.5622
0.0395
SI-IFT-E
0.05830.0007 0.1833
0.0009 0.2671
0.0013 0.4624
0.0022 1.0590
0.0041 2.6394
0.0144
Method15

PMaF: Deep Declarative Layers for Principal Matrix Features

Table 10. Implicit eigen decomposition (IED) layer, GPU memory (MB) for nonsymmetric A. “AutoDiff”: PyTorch eigh() function;“DDN”: deep declarative network; “IFT”: implicit function theorem; “unroll”: unrolling the forward iteration for gradients via PyTorch
autodiff mechanism; “J”: autodiff Jacobian without exploited structure; “E”: our implicit differentiation with exploited structure. Our
best suggestions are highlighted considering the overall solution precision and computational requirements.
5×32
5×64
5×128
5×256
5×512
5×1024
Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward Forward Backward
AutoDiff-unroll
0.5747 0.0972 0.8091 0.3896 1.7515 1.5605 5.5142 6.2456 20.5396 24.9907 80.5903
99.9810
AutoDiff-DDN-E 0.5747 0.0884 0.8076 0.3467 1.7490 1.3789 5.5093 5.5063 20.5298 22.0112 80.5708
88.0210
PI-unroll
0.8071 0.0537 0.9653 0.2280 1.3989 0.9282 2.8350 3.7334
7.5820 14.9688 24.5762
59.9395
PI-DDN-J
0.5137 0.6743 0.5742 5.1694 0.8125 40.6440 1.7603 322.5327
5.5308 2570.0620 20.5718 20520.1206
PI-IFT-J
0.5137 1.2705 0.5742 10.0791 0.8125 80.3135 1.7603 641.2510
5.5308 5125.0010
PI-DDN-E
0.5137 0.0674 0.5742 0.2534 0.8125 1.0044 1.7603 4.0068
5.5308 16.0117 20.5718
64.0215
PI-IFT-E
0.5137 0.1230 0.5742 0.4868 0.8125 1.9409 1.7603 7.7559
5.5308 31.0107 20.5718 124.0205
SI-unroll
4.6055 0.0962 16.5713 0.3882 64.9946 1.5576 256.4995 6.2402 1022.0093 24.9805 4083.0288
99.9609
SI-DDN-J
0.7383 0.6743 1.1025 5.1694 3.1196 40.6440 8.9995 322.5327 32.0093 2570.0620 123.0288 20520.1206
SI-IFT-J
0.7383 1.2705 1.1025 10.0791 3.1196 80.3135 8.9995 641.2510 32.0093 5125.0010
SI-DDN-E
0.7383 0.0674 1.1025 0.2534 3.1196 1.0044 8.9995 4.0068 32.0093 16.0117 123.0288
64.0215
SI-IFT-E
0.7383 0.1230 1.1025 0.4868 3.1196 1.9409 8.9995 7.7559 32.0093 31.0107 123.0288 124.0205
Method16