11institutetext: Interactive Optimization and Learning, Zuse Institute Berlin,
Germany 22institutetext: HTW Berlin, Germany
22email: {mexi, besancon, bolusani, chmiela, hoen<EMAIL_ADDRESS>
# Scylla: a matrix-free fix-propagate-and-project heuristic for mixed-integer
optimization

Gioni Mexi∗ 11 0000-0003-0964-9802 Mathieu Besançon 11 0000-0002-6284-3033
Suresh Bolusani 11 0000-0002-5735-3443 Antonia Chmiela 11 0000-0002-4809-2958
Alexander Hoen 11 0000-0003-1065-1651 Ambros Gleixner 1122 0000-0003-0391-5903

###### Abstract

We introduce Scylla, a primal heuristic for mixed-integer optimization
problems. It exploits approximate solves of the Linear Programming relaxations
through the matrix-free Primal-Dual Hybrid Gradient algorithm with specialized
termination criteria, and derives integer-feasible solutions via fix-and-
propagate procedures and feasibility-pump-like updates to the objective
function. Computational experiments show that the method is particularly
suited to instances with hard linear relaxations.

###### Keywords:

Mixed-integer optimization heuristics matrix-free.

##  1 Introduction

11footnotetext: Corresponding author.

Mixed-integer optimization (MIP) is an essential and challenging class of
problems found in numerous applications including transport, power systems,
engineering, or manufacturing. Primal heuristics are used as stand-alone
methods to compute good-quality solutions or as subroutines in branch-and-
bound solvers, potentially exploiting the information extracted in the rest of
the solution process [5]. Solutions derived from primal heuristics can be
returned immediately when the corresponding gap is small, they help to prune
nodes of the branch-and-bound tree, and are also leveraged by other
subroutines of the solver, e.g., cut selection [12].

When integrated within a branch-and-bound scheme, many primal heuristics
assume that an optimal solution of the linear programming (LP) relaxation is
readily available and that the LP can be re-optimized quickly after few
modifications. These repeated calls to LP solvers can become a bottleneck in
several scenarios. First, while it is true that standard MIP solvers produce
and rely on optimal LP relaxation solutions, they are not yet available during
presolving. Second, for some applications, re-optimization via the simplex
method is expensive [14]. For primal heuristics in a stand-alone setting, the
computation of an initial solution to the LP relaxation by standard methods
that require matrix factorizations, such as the simplex or interior point
methods, may be expensive for large-scale applications.

To address this problem, we present _Scylla_ , a _matrix-free_ primal
heuristic, in the sense that it uses only iterated matrix-vector
multiplications to compute approximate LP solutions without factorizing the
constraint matrix nor directly solving linear systems. We supplement this with
LP-free constraint propagation to round such approximate reference solutions
to integer feasibility and objective function updates inspired by the well-
known feasibility pump [9, 4, 2].

The paper is organized as follows. In Section 2, we present the algorithm
design and highlight its key features. In Section 3, we report on our
computational experiments. We summarize our observations and discuss future
research directions in Section 4.

##  2 The algorithm

The overall structure of Scylla is outlined in Algorithm 0.1 and is based upon
several subroutines:

  * •

a limited-precision call to the Primal-Dual Hybrid Gradient (PDHG) algorithm
[17, 8] to compute an approximate solution of the LP relaxation,

  * •

a fix-and-propagate procedure to attempt rounding this fractional solution to
a fully, i.e., also integer, feasible solution,

  * •

an objective function update similar to that of the feasibility pump [9, 4, 2]
in order to iteratively drive the PDHG solution towards integer feasibility.

The outer loop of the algorithm is relatable to a feasibility pump [9, 4]
where the rounding step uses constraint propagation [10] and the projection
step is replaced with a limited number of PDHG iterations until relaxed
termination criteria are met. The inner iterations counted by $k$ are
dominated by matrix-vector multiplications with the constraint matrix and its
transpose; the outer iterations counted by $K$ are dominated by a fix-and-
propagate call and update the objective function at the end. Note that before
starting the actual algorithm, we presolve the instance to reduce its size
using the SCIP solver [7].

In the following sections, we describe the principal components of the
algorithm and their interaction.

Input : MIP instance $\min\\{c^{\top}x\,:\,Ax\geq b,\ \ell\leq x\leq u,\
x_{j}\in\mathbb{Z},\ \forall\ j\in I\\}$

Output : MIP-feasible solution $x$ or $\varnothing$ if none found

1 initialize $k=1$, $K=1$, $\alpha=0.9$

2 initialize primal-dual step sizes $\tau_{x}>0,\tau_{y}>0$ and
$\bar{x}^{(k)}=0$, $\bar{y}^{(k)}=0$

3 initialize $\hat{c}=c$

4 while  _limits not reached_ do

5 repeat

update step sizes $\tau_{x},\tau_{y}$

/* PDHG iterations */

6 $d_{x}=A^{\top}\bar{y}^{(k)}-\hat{c}$ and
$\bar{x}^{(k+1)}=\texttt{proj}_{[\ell,u]}(\bar{x}^{(k)}+\tau_{x}d_{x})$

7 $d_{y}=b-A(2\bar{x}^{(k+1)}-\bar{x}^{(k)})$ and
$\bar{y}^{(k+1)}=\texttt{proj}_{\mathbb{R}^{n}_{+}}(\bar{y}^{(k)}+\tau_{y}d_{y})$

8 $k=k+1$

9

10 until  _$\bar{x}^{(k)},\bar{y}^{(k)}\approx$ LP optimal or bounds stalling
or $k>\mathtt{iterlimit}$_

11 if  _$\bar{x}^{(k)}$ MIP-feasible_ then return _$\bar{x}^{(k)}$_

$\hat{x}^{(K)}=\texttt{fix-and-propagate}(\bar{x}^{(k)})$

/* rounding attempt */

12 if  _$\hat{x}^{(K)}$ MIP-feasible_ then return _$\hat{x}^{(K)}$_

13 if  _cycling on $\hat{x}^{(K)}$ detected_ then
$\hat{x}^{(K)}=\texttt{perturb}(\hat{x}^{(K)}$)

$\hat{c}^{\top}x=\alpha^{K}\frac{\sqrt{\lvert I\rvert}}{\left\lVert
c\right\rVert}c^{\top}x+(1-\alpha^{K})\Delta(x,\hat{x}^{(K)})$

/* objective update */

14 $K=K+1$

return _$\varnothing$_

Algorithm 0.1 High-level structure of the Scylla heuristic

###  2.1 Approximate LP solutions by matrix-free linear optimization

Heuristic and exact methods for mixed-integer problems often rely on repeated
calls to a linear optimization solver. These calls represent a substantial
amount of the runtime and require, for the most common LP solving techniques,
a matrix factorization. As highlighted in [3], first-order methods applied to
the saddle-point formulations of LPs open the door to a trade-off between
solution accuracy and runtime which was rarely reasoned about in the linear
and mixed-integer setting. In our algorithm, these approximate LP solutions
are used as a starting point for a fix-and-propagate procedure.

To leverage this trade-off, we use the PDHG [8] algorithm and its recent PDLP
implementation (PDHG specialized for linear optimization) [3]. Some notable
algorithmic improvements were introduced in [3] compared to the standard PDHG
including adaptive step size strategies, restarts, and matrix preconditioning.
We use the PDLP implementation available in OR-Tools111The package is
available at https://github.com/google/or-tools.. The modifications allow us
to warm-start the solution process from a previously computed primal-dual pair
and to add custom stopping criteria detailed in Section 2.2.

For conciseness of the algorithm description, we only included the core PDHG
iterations and did not detail the step size computation procedure or
restarting scheme from PDLP. For Lines 0.1 to 0.1 of Algorithm 0.1, Scylla
calls PDLP to solve similar LPs where only the objective function vector
changes. We use the fact that the state of PDLP is fully described by its
primal-dual pair and primal and dual residuals to warm-start all calls after
the first one, which may help to reduce the total number of iterations after
the first solve. Using PDLP also means that Scylla remains a matrix-free
heuristic only requiring sparse matrix-vector multiplications.

###  2.2 Custom stopping criteria

The standard stopping criterion for PDLP is a maximum error $\varepsilon$ on
the primal and dual feasibilities. This error can be relaxed compared to what
would be expected in “standard” LP solving since the resulting fractional
solution is then guiding a rounding procedure. In particular, this error can
also be adapted dynamically to allow a bigger tolerance at the beginning of
the heuristic. In addition to increasing the initial error tolerance (to
$\varepsilon^{(1)}=0.01$ in our experiments), we also use a progressive
refinement scheme in which the error is reduced geometrically by a factor
$\beta=0.98$, resulting in the following.

$\varepsilon^{(K+1)}=\max\\{\beta\ \varepsilon^{(K)},10^{-8}\\}$

###  2.3 The fix-and-propagate procedure

Fix-and-propagate algorithms have been used successfully as part of branch-
and-bound solvers [6, 11]. They are also an essential step in the most
successful variants of feasibility pump algorithms [10], and we use them with
the same motivation: to “round” the PDHG solution $\bar{x}^{(k)}$ to integer
feasibility.

“Rounding” is achieved by fixing an integer variable with fractional value to
an integer within its domain, followed by a propagating step. This is repeated
until all such variables are fixed or infeasibility is detected by some domain
becoming empty. In the latter case, fix-and-propagate continues in order to
produce an integer vector by ignoring any constraint that would lead to empty
domains. At the end of the propagation, all remaining unfixed variables are
fixed to their values in the fractional reference solution or its projection
on to their domain. The procedure always produces an integer-feasible, but not
necessarily LP-feasible, solution.

One major issue of feasibility pump-like algorithms and hence also Scylla is
cycling: after some iterations, it may happen that
$\hat{x}^{(K)}=\hat{x}^{(K^{\prime})}$ for $1\leq K^{\prime}<K$. Similar to
[9], this issue is handled by randomly perturbing some of the variables in the
current rounded vector $\hat{x}^{(K)}$ (Line 0.1).

###  2.4 Objective function update

Inspired by objective feasibility pump approaches [9, 4], we update the
objective function as follows. After computing the integer point
$\hat{x}^{(K)}$ and before the subsequent PDLP call, the objective function is
updated (Line 0.1 of Algorithm 0.1) as a convex combination of the initial
objective function $c^{\top}x$ and an $\ell_{1}$-norm distance to the integer
vector $\Delta(x,\hat{x}^{(K)})$ with $\alpha\in[0,1)$ chosen as constant
$0.9$ in our experiments.

This convex combination represents a trade-off between the initial objective
(towards high-quality solutions) and the distance to an integer vector
(towards integer feasibility). After each iteration $K$, the weight shifts
further towards integer feasibility in the attempt to produce a PDLP solution
$\bar{x}^{(k)}$ that fix-and-propagate can successfully convert to a MIP-
feasible solution.

##  3 Computational results

As a baseline for our implementation, we modified an open source
implementation of the Feasibility Pump [16]
(FP)222https://github.com/GioniMexi/FeasPumpCollection (that uses SoPlex 6.0.2
[7] as the underlying LP solver) and extended it to support approximate LP
solves with PDLP.

Our test set consists of the benchmark library MIPLIB 2017 [13] and the
instances used in [10]. We excluded all instances that were solved to
optimality during presolving with SCIP, leaving us with 283 instances.

All results were obtained on a cluster of Intel Xeon Gold 6338 @ 2.0 GHz
machines with 1024 GB of RAM and a time limit of $3600$ seconds. To mitigate
the impact of performance variability [15], both FP and Scylla are tested on
the complete test set with ten different random seeds.

Table 1 shows the summarized results for 2830 instances (283 instances
$\times$ 10 random seeds). On the whole testset, FP outperforms Scylla w.r.t.
both solution quality and running time. This is an expected result since the
LPs in the FP iterations are for most instances easily solved to optimality,
and warm-starting of simplex LP solvers is very effective. Scylla therefore
loses the benefit of exact LP solutions while still requiring similar costs
per iteration. We see this when looking at the two different subsets of
instances.

For the instances where the root LP solves in $[1,60]$ seconds, FP is able to
find more solutions than Scylla and provides better quality solutions. Scylla
is however faster on average, despite suffering under a greater model
construction overhead which increases with the instance size. Note that this
overhead is tied to the current implementation and not to the heuristic
itself. More specifically, the current OR-Tools interface through which PDLP
is accessed requires re-building the model at every iteration of the
feasibility pump, a restriction that penalized our approach, in particular for
larger instances.

For instances in $[60,\text{lim}]$, Scylla significantly improves over FP in
terms of both the number of solved instances and running time. Notably, Scylla
is capable of finding feasible solutions for 42 more instances than FP, and it
is nearly four times faster. These results suggest that for MIPs with
difficult LP relaxations, Scylla is the preferred heuristic choice to compute
first primal solutions.

Table 1: Aggregated results comparing FP with Scylla. Row $[1,60]$ ($[60,\text{lim}]$) only considers instances on which solving the LP takes between 1 and 60 seconds (more than 60 seconds). Column ‘#’ (‘only’) refers to the number of instances for which a solution was found (exclusively by only one heuristic), ‘gap’ refers to the relative gap to the best known solution, ‘wins’ refers to the number of instances with a performance improvement of 10% over the other heuristic, and ‘overhead’ refers to the overhead time to build the model. |  | found sol | #wins | shift. geom. mean |
---|---|---|---|---|---
testset | algorithm | # | only | gap | time(s) | gap | time(s) | overhead(s)
all | FP | 2268 | 152 | 1154 | 1557 | 27.46 | 10.06 | 0.02
| Scylla | 2270 | 154 | 485 | 797 | 40.92 | 17.30 | 1.93
$[1,60]$ | FP | 660 | 62 | 341 | 295 | 26.75 | 43.16 | 0.01
| Scylla | 635 | 37 | 130 | 381 | 41.78 | 36.75 | 1.99
$[60,\text{lim}]$ | FP | 163 | 40 | 98 | 64 | 28.11 | 974.66 | 0.01
| Scylla | 205 | 82 | 105 | 180 | 35.37 | 234.80 | 6.09

##  4 Conclusion

In this paper, we presented Scylla, a primal heuristic for mixed-integer
linear optimization that combines matrix-free computation of approximate LP
solutions with a fix-and-propagate procedure and an objective function update
scheme similar to the feasibility pump. Computational experiments on a large
set of instances from the MIPLIB 2017 collection show that the method is
particularly suited to instances with hard linear relaxations. Our current
implementation does not exploit the potential for parallelization that lies in
the sparse matrix-vector multiplications of PDHG, which dominate the PDLP
runtime. A robust, efficient, and multithreaded sparse BLAS library will thus
accelerate this core operation of PDLP and could be key to a broader success
especially on large-scale instances.

One aspect to be explored in future research is the fact that PDHG is by
design suited to solve more general convex differentiable optimization
problems, and PDLP handles quadratic problems where the objective Hessian is
diagonal [1]. On the one hand, this makes Scylla a candidate heuristic for the
broader class of MIQPs. On the other hand, this can be exploited also in the
linear setting by choosing the $\ell_{2}$-norm as the distance function
$\Delta(\cdot,\hat{x}^{(K)})$ and avoiding auxiliary variables needed to model
the $\ell_{1}$-norm on general integer variables [9, 4].

## Acknowledgements

Research reported in this paper was partially supported through the Research
Campus Modal funded by the German Federal Ministry of Education and Research
(fund numbers 05M14ZAM,05M20ZBM) and the Deutsche Forschungsgemeinschaft (DFG)
through the DFG Cluster of Excellence MATH+.

## References

  * [1] Mathematical background for PDLP, Google OR-tools documentation. https://developers.google.com/optimization/lp/pdlp_math, accessed: 2023-06-27
  * [2] Achterberg, T., Berthold, T.: Improving the feasibility pump. Discrete Optimization 4(1), 77–86 (2007)
  * [3] Applegate, D., Díaz, M., Hinder, O., Lu, H., Lubin, M., O’Donoghue, B., Schudy, W.: Practical large-scale linear programming using primal-dual hybrid gradient. Advances in Neural Information Processing Systems 34 (2021)
  * [4] Bertacco, L., Fischetti, M., Lodi, A.: A feasibility pump heuristic for general mixed-integer problems. Discrete Optimization 4(1), 63–76 (2007)
  * [5] Berthold, T.: Heuristic algorithms in global MINLP solvers. Ph.D. thesis, Technische Universität Berlin (2014)
  * [6] Berthold, T., Hendel, G.: Shift-and-propagate. Journal of Heuristics 21(1), 73 – 106 (2014). https://doi.org/10.1007/s10732-014-9271-0
  * [7] Bestuzheva, K., Besançon, M., Chen, W.K., Chmiela, A., Donkiewicz, T., van Doornmalen, J., Eifler, L., Gaul, O., Gamrath, G., Gleixner, A., et al.: Enabling research through the scip optimization suite 8.0. ACM Transactions on Mathematical Software (2023)
  * [8] Chambolle, A., Pock, T.: A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of mathematical imaging and vision 40(1), 120–145 (2011)
  * [9] Fischetti, M., Glover, F., Lodi, A.: The feasibility pump. Mathematical Programming 104(1), 91–104 (2005)
  * [10] Fischetti, M., Salvagnin, D.: Feasibility pump 2.0. Mathematical Programming Computation 1(2), 201–222 (2009)
  * [11] Gamrath, G., Berthold, T., Heinz, S., Winkler, M.: Structure-driven fix-and-propagate heuristics for mixed integer programming. Mathematical Programming Computation 11(4), 675–702 (2019)
  * [12] Gleixner, A., Bastubbe, M., Eifler, L., Gally, T., Gamrath, G., Gottwald, R.L., Hendel, G., Hojny, C., Koch, T., Lübbecke, M.E., Maher, S.J., Miltenberger, M., Müller, B., Pfetsch, M.E., Puchert, C., Rehfeldt, D., Schlösser, F., Schubert, C., Serrano, F., Shinano, Y., Viernickel, J.M., Walter, M., Wegscheider, F., Witt, J.T., Witzig, J.: The SCIP Optimization Suite 6.0. Tech. rep., Optimization Online (2018)
  * [13] Gleixner, A., Hendel, G., Gamrath, G., Achterberg, T., Bastubbe, M., Berthold, T., Christophel, P., Jarck, K., Koch, T., Linderoth, J., et al.: MIPLIB 2017: data-driven compilation of the 6th mixed-integer programming library. Mathematical Programming Computation pp. 1–48 (2021). https://doi.org/10.1007/s12532-020-00194-3
  * [14] Koch, T., Achterberg, T., Andersen, E., Bastert, O., Berthold, T., Bixby, R.E., Danna, E., Gamrath, G., Gleixner, A.M., Heinz, S., Lodi, A., Mittelmann, H., Ralphs, T., Salvagnin, D., Steffy, D.E., Wolter, K.: MIPLIB 2010. Mathematical Programming Computation 3(2), 103–163 (2011). https://doi.org/10.1007/s12532-011-0025-9
  * [15] Lodi, A., Tramontani, A.: Performance variability in mixed-integer programming. In: Theory Driven by Influential Applications, pp. 1–12. INFORMS (2013). https://doi.org/10.1287/educ.2013.0112
  * [16] Mexi, G., Berthold, T., Salvagnin, D.: Using multiple reference vectors and objective scaling in the feasibility pump. EURO Journal on Computational Optimization 11, 100066 (2023)
  * [17] Zhu, M., Chan, T.: An efficient primal-dual hybrid gradient algorithm for total variation image restoration. UCLA Cam Report 34, 8–34 (2008)


# Invariant subspaces of the shift plus complex Volterra operator

Željko Čučković, Bhupendra Paudyal  This paper is based on a research which
forms a part of the second author’s University of Toledo Ph.D. Thesis written
under the supervision of Professor Željko Čučković

###### Abstract

The lattice of closed invariant subspaces of the Volterra operator acting on
$L^{2}(0,1)$ was completely described by Sarason [8]. On the other hand, he
explicitly found the lattice of closed invariant subspaces of the shift plus
Volterra operator on $L^{2}(0,1)$ in [9]. Inspired by Sarason’s results, we
find the lattice of closed invariant subspaces of the shift plus complex
Volterra operator acting on the Hardy space.

##  1 Introduction

The invariant subspace problem is a well known problem in operator theory. The
conjecture is that every bounded linear operator $T$ on a separable Hilbert
space $H$ has a non-trivial closed invariant subspace (a closed linear
subspace W of $H$ which is different from $\\{0\\}$ and $H$ such that
$T(W)\subset W$). Over the years it has been shown that many familiar classes
of operators do have invariant subspaces. The goal of this paper is to
characterize the lattice of closed invariant subspaces of the shift plus
complex Volterra operator on the Hardy space.

Let $\mathbb{D}$ be the unit disk of the complex plane and $H(\mathbb{D})$ be
the space of holomorphic functions on the unit disk. We say that a holomorphic
function $f(z)=\sum_{n=0}^{\infty}a_{n}z^{n}$ on the unit disk belongs to the
Hardy space $\mathcal{H}^{2}(\mathbb{D})$, if its sequence of power series
coefficients is square-summable:

$\mathcal{H}^{2}(\mathbb{D})=\\{f\in
H(\mathbb{D}):\sum_{n=0}^{\infty}|a_{n}|^{2}<\infty\\}.$

We define a norm on $\mathcal{H}^{2}(\mathbb{D})$ by

$\|f\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}=\sum_{n=0}^{\infty}|a_{n}|^{2}.$ (1)

It is well known that $\mathcal{H}^{2}(\mathbb{D})$ is a Hilbert space with
the inner product

$\langle
f,g\rangle_{\mathcal{H}^{2}(\mathbb{D})}=\sum_{n=0}^{\infty}a_{n}\overline{b_{n}}$

for $f(z)=\sum_{n=0}^{\infty}a_{n}z^{n}$ and
$g(z)=\sum_{n=0}^{\infty}b_{n}z^{n}$.

The operator defined on $\mathcal{H}^{2}(\mathbb{D})$

$(M_{z}f)(z)=zf(z)\hskip
14.22636pt\text{for}~{}f\in\mathcal{H}^{2}(\mathbb{D})~{}\text{and}~{}z\in\mathbb{D}$

is called the shift operator. The lattice of the shift operator acting on the
Hardy space is completely described by Beurling’s Theorem [2], and it is one
of the the most celebrated and widely used results. Let $L^{2}(0,1)$ be the
space of square integrable functions on $(0,1)$. Sarason [8] characterized all
closed invariant subspaces of the Volterra operator

$(Vf)(x)=\int_{0}^{x}f(y)dy\hskip 14.22636pt\text{for}~{}f\in
L^{2}(0,1)~{}\text{and}~{}0<x\leq y<1.$

Aleman and Korenblum studied the complex Volterra operator in
$\mathcal{H}^{2}(\mathbb{D})$ defined by

$(\mathcal{V}f)(z)=\int_{0}^{z}f(w)dw.$

Then they characterized the lattice of closed invariant subspaces of
$\mathcal{V}$ in [1]. While doing so they used the Beurling’s Theorem. On the
other hand, Sarason [9] studied the lattice of closed invariant subspaces of
multiplication by $x$ plus Volterra operator, $M_{x}+V$ acting on
$L^{2}(0,1)$. Montes-Rodriguez, Ponce-Escudero and Shkarin [6] and Cowen,
Gunatillake and Ko [3] used the idea of Sarason to study the invariant
subspaces of certain classes of composition operators on Hardy spaces.

Following Sarason’s work we are interested in characterizing the lattice of
closed invariant subspaces of the shift plus complex Volterra operator on the
Hardy space. Denote by $T$ the operator

$(Tf)(z)=zf(z)+\int_{0}^{z}f(w)dw,\hskip
14.22636pt\text{for}~{}f\in\mathcal{H}^{2}(\mathbb{D})~{}\text{and}~{}z\in\mathbb{D}.$
(2)

Since the shift operator is an isometry and the complex Volterra operator is a
contraction, $T$ is clearly bounded operator on the Hardy space. To prove our
main result we use the space $\mathcal{S}^{2}(\mathbb{D})$ defined by

$\mathcal{S}^{2}(\mathbb{D})=\\{f\in
H(\mathbb{D}):Df\in\mathcal{H}^{2}(\mathbb{D})\\},$

where $D=\frac{d}{dz}$ is the differential operator. It is clear that if $Df$
is in $\mathcal{H}^{2}(\mathbb{D})$, then $f$ belongs to
$\mathcal{H}^{2}(\mathbb{D})$. The norm of $\mathcal{S}^{2}(\mathbb{D})$ is
defined by

$\|f\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}=\|Df\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+\|f\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}.$
(3)

Corresponding inner product is given by

$\langle f,g\rangle_{\mathcal{S}^{2}(\mathbb{D})}=\langle
Df,Dg\rangle_{\mathcal{H}^{2}(\mathbb{D})}+\langle
f,g\rangle_{\mathcal{H}^{2}(\mathbb{D})}.$

Our work describes the lattice of closed invariant subspaces of $T$ acting on
$\mathcal{H}^{2}(\mathbb{D})$. The main tool in our proof is Korenbljum’s
result from [5], in which he characterized all closed ideals of
$\mathcal{S}^{2}(\mathbb{D})$. Korenbljum’s paper, in turn, relates to the
work of Rudin [7] and his main theorem on the structure of the closed ideals
of the disk algebra. The following result is known to the specialists but we
could not find the proof in the literature. We include the proof for the sake
of completeness.

######  Proposition 1.

The following statements are true:

  1. i.

$\mathcal{S}^{2}(\mathbb{D})\subset H^{\infty}.$

  2. ii.

$\mathcal{S}^{2}(\mathbb{D})$ is Banach algebra.

  3. iii.

Polynomials are dense in $\mathcal{S}^{2}(\mathbb{D})$.

###### Proof.

  1. i.

Let $f\in\mathcal{S}^{2}(\mathbb{D})$, then $Df\in\mathcal{H}^{2}(\mathbb{D})$
and hence

$|Df(z)|\leq\dfrac{\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}}{\sqrt{1-|z|^{2}}},~{}\hskip
14.22636pt~{}\text{for~{} all}~{}z\in\mathbb{D}.$

Now,

$\displaystyle|f(z)-f(0)|=$ $\displaystyle|\int_{0}^{1}zDf(tz)dt|$
$\displaystyle\leq$ $\displaystyle\int_{0}^{1}|zDf(tz)|dt$ $\displaystyle\leq$
$\displaystyle\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}\int_{0}^{1}\frac{|z|}{\sqrt{1-|tz|^{2}}}dt$
$\displaystyle\leq$
$\displaystyle\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}\int_{0}^{1}\frac{|z|}{\sqrt{(1-|tz|)(1+t|z|)}}dt.$

Since $1+t|z|\geq 1$, we have $\dfrac{1}{\sqrt{1+t|z|}}\leq 1$.

$\displaystyle|f(z)-f(0)|\leq$
$\displaystyle\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}\int_{0}^{1}\frac{|z|}{\sqrt{1-|tz|}}dt$
$\displaystyle\leq$ $\displaystyle 2\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}.$ (4)

This clearly shows that $f$ belongs to $H^{\infty}$, and hence
$\mathcal{S}^{2}(\mathbb{D})\subset H^{\infty}.$

  2. ii.

To show $\mathcal{S}^{2}(\mathbb{D})$ is a Banach space under the above norm,
suppose $\\{g_{n}\\}_{0}^{\infty}$ is Cauchy in $\mathcal{S}^{2}(\mathbb{D})$
norm. Then, clearly $\\{g_{n}\\}$ and $\\{Dg_{n}\\}$ are Cauchy on
$\mathcal{H}^{2}(\mathbb{D})$ norm. Since $\mathcal{H}^{2}(\mathbb{D})$ is a
Banach space, $\\{g_{n}\\}$ converges to a holomorphic function
$g\in\mathcal{H}^{2}(\mathbb{D})$ and $\\{Dg_{n}\\}$ converges to the function
$Dg\in\mathcal{H}^{2}(\mathbb{D})$. Therefore
$g\in\mathcal{S}^{2}(\mathbb{D})$, and hence $\mathcal{S}^{2}(\mathbb{D})$ is
a Banach space under the given norm. Pointwise multiplication on
$\mathcal{S}^{2}(\mathbb{D})$ form an algebra. For this, suppose that $f$ and
$g$ are in $\mathcal{S}^{2}(\mathbb{D})$.

$\displaystyle\|fg\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}=$
$\displaystyle\|D(fg)\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+\|fg\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}$
$\displaystyle=$
$\displaystyle\|gDf+fDg\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+\|fg\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}$
$\displaystyle\leq$
$\displaystyle\|gDf\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+2\|gDf\|_{\mathcal{H}^{2}(\mathbb{D})}\|fDg\|_{\mathcal{H}^{2}(\mathbb{D})}+\|fDg\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+\|fg\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}.$
(5)

Using (i), we see that for any $f\in\mathcal{S}^{2}(\mathbb{D})$

$\displaystyle\|f\|_{\infty}\leq$ $\displaystyle
2\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}+|f(0)|$ $\displaystyle\leq$
$\displaystyle
2\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}+2\|f\|_{\mathcal{H}^{2}(\mathbb{D})}=2\|f\|_{\mathcal{S}^{2}(\mathbb{D})}.$

Hence using (ii),

$\displaystyle\|fg\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}\leq$
$\displaystyle\|g\|^{2}_{\infty}\|Df\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+2\|f\|_{\infty}\|g\|_{\infty}\|Df\|_{\mathcal{H}^{2}(\mathbb{D})}\|Dg\|_{\mathcal{H}^{2}(\mathbb{D})}$
$\displaystyle+\|f\|^{2}_{\infty}\|Dg\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+\|g\|^{2}_{\infty}\|f\|^{2}_{\infty}$
$\displaystyle\leq$ $\displaystyle
4\|g\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}\|f\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}+8\|g\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}\|f\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}+4\|g\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}\|f\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}$
$\displaystyle\leq$ $\displaystyle
16\|f\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}\|g\|^{2}_{\mathcal{S}^{2}(\mathbb{D})}.$

######  Remark 1.

Clearly $\mathcal{S}^{2}(\mathbb{D})$ is an algebra with unit and with a norm
$\|.\|_{\mathcal{S}^{2}(\mathbb{D})}$ under which it is a Banach space.
Furthermore, we see that the multiplication is continuous in each factor
separately. Therefore there exists a norm equivalent to
$\|.\|_{\mathcal{S}^{2}(\mathbb{D})}$, for which $\mathcal{S}^{2}(\mathbb{D})$
is a Banach algebra (see [4, page 212]). Keeping in mind the fact above, we do
not differentiate the norm $\|.\|_{\mathcal{S}^{2}(\mathbb{D})}$ and its
equivalent norm that makes $\mathcal{S}^{2}(\mathbb{D})$ a Banach algebra.

  3. iii.

We want to show polynomials are dense in $\mathcal{S}^{2}(\mathbb{D})$. Given
$f\in\mathcal{S}^{2}(\mathbb{D})$, let $f_{q}(z)=f(qz)$ be its dilation and
$Df_{q}(z)=qDf(qz)=q(Df)_{q}(z)$ be derivative of dilation where $0<q<1$. Each
function $f_{q}$ is analytic in a larger disk, so it can be approximated
uniformly on $\mathbb{D}$ by a sequence of holomorphic polynomials
$P_{q}^{n}$, and hence $Df_{q}$ can be approximated uniformly on $\mathbb{D}$
by holomorphic polynomials $DP_{q}^{n}$. So it will be enough to prove that
$f$ can be approximated in $\mathcal{S}^{2}(\mathbb{D})$ by its dilation. That
is to say $\|f-f_{q}\|_{\mathcal{S}^{2}(\mathbb{D})}\rightarrow 0$ as
$q\rightarrow 1$. This means that
$\|f-f_{q}\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}+\|Df-
Df_{q}\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}\rightarrow 0$ as $q\rightarrow 1$.
So finally it is enough to prove that
$\|f-f_{q}\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}\rightarrow 0$ and $\|Df-
Df_{q}\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}\rightarrow 0$ as $q\rightarrow 1$.
For this, let us assume

$\displaystyle{f(z)=\sum_{n=0}^{\infty}a_{n}z^{n}}.$

Since $f\in\mathcal{H}^{2}(\mathbb{D})$, for all $\epsilon>0$, we can choose a
natural number $N$ large enough such that

$\displaystyle{\sum_{n=N+1}^{\infty}|a_{n}|^{2}<\frac{\epsilon}{2}}.$

Now choose $q_{\epsilon}\in(0,1)$ such that

$(1-q_{\epsilon}^{N})^{2}\sum_{n=0}^{N}|a_{n}|^{2}<\frac{\epsilon}{2}.$

Then, since

$\|f-f_{q}\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}=\left\|\sum_{n=0}^{\infty}a_{n}z^{n}(1-q^{n})\right\|^{2}=\displaystyle{\sum_{n=0}^{\infty}|a_{n}(1-q^{n})|^{2}},$

it follows that for all $q\geq q_{\epsilon}$

$\displaystyle\|f-f_{q}\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}==$
$\displaystyle\sum_{n=0}^{N}|a_{n}(1-q^{n})|^{2}+\sum_{n=N+1}^{\infty}|a_{n}(1-q^{n})|^{2}$
$\displaystyle\leq$
$\displaystyle(1-q^{N})\sum_{n=0}^{N}|a_{n}|^{2}+\sum_{n=N+1}^{\infty}|a_{n}|^{2}$
$\displaystyle\leq$
$\displaystyle\frac{\epsilon}{2}+\frac{\epsilon}{2}=\epsilon.$

This shows that $\|f-f_{q}\|^{2}_{\mathcal{H}^{2}(\mathbb{D})}\rightarrow 0$
as $q\rightarrow 1$. On the other hand, we have

$\displaystyle|Df(z)-Df_{q}(z)|=$
$\displaystyle\left|\sum_{n=1}^{\infty}na_{n}z^{n-1}-\sum_{n=1}^{\infty}na_{n}q^{n}z^{n-1}\right|=$
$\displaystyle\left|\sum_{n=1}^{\infty}na_{n}(1-q^{n})z^{n-1}\right|.$

Similarly, we can show that $\|Df-Df_{q}\|^{2}\rightarrow 0$ as $q$ approaches
to $1$.

∎

######  Definition 1.

Define

$\mathcal{S}_{0}^{2}(\mathbb{D})=\\{f\in\mathcal{S}^{2}(\mathbb{D}):f(0)=0\\}$

######  Corollary 1.

$\mathcal{S}_{0}^{2}(\mathbb{D})\subset\mathcal{S}^{2}(\mathbb{D})$ is a
Banach algebra with the norm defined for $\mathcal{S}^{2}(\mathbb{D})$ and
$\mathcal{S}^{2}(\mathbb{D})=[1]\bigoplus\mathcal{S}_{0}^{2}(\mathbb{D})$, and
hence

$\mathcal{S}_{0}^{2}(\mathbb{D})=\overline{span}\\{z^{n}:n\in\mathbb{N}\\}.$

###### Proof.

For any $f$ and $g$ in
$\mathcal{S}_{0}^{2}(\mathbb{D})\subset\mathcal{S}^{2}(\mathbb{D})$, we
immediately see that

$\|fg\|_{\mathcal{S}^{2}(\mathbb{D})}\leq\|f\|_{\mathcal{S}^{2}(\mathbb{D})}\|g\|_{\mathcal{S}^{2}(\mathbb{D})}.$

Also, we have $f(0)=0$ and $g(0)=0$ so $(fg)(0)=0$ and hence $fg$ belongs to
$\mathcal{S}_{0}^{2}(\mathbb{D})$. To show $\mathcal{S}_{0}^{2}(\mathbb{D})$
is a closed subalgebra of $\mathcal{S}^{2}(\mathbb{D})$, assume
$g\in\overline{\mathcal{S}_{0}^{2}(\mathbb{D})}$. That means there exists a
sequence $g_{n}\in\mathcal{S}_{0}^{2}(\mathbb{D})$, $n\in\mathbb{N}$ such that
$g_{n}$ converges to $g$ in $\mathcal{S}^{2}(\mathbb{D})$ norm. This implies
$g_{n}$ converges to $g$ in $\mathcal{H}^{2}(\mathbb{D})$ norm. Since
$g_{n}(0)=0$ for all $n$, it follows that $g(0)=0$. Moreover, the other part
of the corollary follows from the Definition 1. ∎

The following theorem plays important role in the proof of main theorem.

######  Theorem 1.

Let $\mathcal{V}$ be the Volterra operator on $\mathcal{H}^{2}(\mathbb{D})$.
Then the following statements are true:

  1. i.

Range of $\mathcal{V}=\mathcal{S}_{0}^{2}(\mathbb{D})$.

  2. ii.

$\mathcal{V}$ is a bounded isomorphism from $\mathcal{H}^{2}(\mathbb{D})$ onto
$\mathcal{S}_{0}^{2}(\mathbb{D})$, and its inverse is $D$.

  3. iii.

The operator $T$ acting on $\mathcal{H}^{2}(\mathbb{D})$ is similar under $V$
to the multiplication operator $M_{z}$ acting on
$\mathcal{S}_{0}^{2}(\mathbb{D})$.

###### Proof.

  1. i.

Let $g$ be in the range of $\mathcal{V}$, then there exists
$f\in\mathcal{H}^{2}(\mathbb{D})$ such that

$\displaystyle{g(z)=(\mathcal{V}f)(z)=\int_{0}^{z}f(w)dw}$

then $Dg=f\in\mathcal{H}^{2}(\mathbb{D})$ and $g(0)=0$. Hence
$g\in\mathcal{S}_{0}^{2}(\mathbb{D})$. Conversely, suppose that $g$ belongs to
$\mathcal{S}_{0}^{2}(\mathbb{D})$, then

$\displaystyle{(\mathcal{V}Dg)(z)=\int_{0}^{z}(Dg)(w)dw=g(z)-g(0)=g(z)}.$

Therefore $g$ belongs to the range of $\mathcal{V}$.

  2. ii.

First we want to show $\mathcal{V}$ is a bounded operator on
$\mathcal{H}^{2}(\mathbb{D})$. Let us assume $f$ is in the Hardy space.

$\displaystyle\|\mathcal{V}(f)\|_{\mathcal{S}^{2}(\mathbb{D})}=$
$\displaystyle\|\mathcal{V}(f)\|_{\mathcal{H}^{2}(\mathbb{D})}+\|D(\mathcal{V}(f))\|_{\mathcal{H}^{2}(\mathbb{D})}$
$\displaystyle\leq$
$\displaystyle\|f\|_{\mathcal{H}^{2}(\mathbb{D})}+\|f\|_{\mathcal{H}^{2}(\mathbb{D})}$
$\displaystyle=$ $\displaystyle 2\|f\|_{\mathcal{H}^{2}(\mathbb{D})}.$

Hence the map $\mathcal{V}$ from $\mathcal{H}^{2}(\mathbb{D})$ onto
$\mathcal{S}_{0}^{2}(\mathbb{D})$ is bounded. Clearly $\mathcal{V}$ is linear.
Now to show $\mathcal{V}$ is one-one, assume that $f_{1}$ and $f_{2}$ belong
the Hardy space, and also assume that

$\int_{0}^{z}f_{1}(w)dw=\int_{0}^{z}f_{2}(w)dw,\hskip
28.45274pt\text{for~{}all}~{}z\in\mathbb{D}.$

Differentiating both sides we see that $f_{1}=f_{2}$ and hence $\mathcal{V}$
is one-one. From part $(i)$ we have $\mathcal{V}Dg=g$ and clearly
$D\mathcal{V}f=f$. This shows that $\mathcal{V}$ is a bounded bijective linear
operator from $\mathcal{H}^{2}(\mathbb{D})$ onto
$\mathcal{S}_{0}^{2}(\mathbb{D})$ and $\mathcal{V}^{-1}=D$.

  3. iii.

Suppose $f$ belongs to $\mathcal{H}^{2}(\mathbb{D})$ and also suppose
$\mathcal{V}f=g$, for some $g\in\mathcal{S}_{0}^{2}(\mathbb{D})$. Therefore we
have $f(z)=(\mathcal{V}^{-1}g)(z)=(Dg)(z)$.

$\displaystyle(Tf)(z)=$ $\displaystyle zf(z)+(\mathcal{V}f)(z)$
$\displaystyle=$ $\displaystyle z(Dg)(z)+g(z)$ $\displaystyle=$ $\displaystyle
D(zg(z)).$

Now applying $\mathcal{V}$ on the both side, we see that

$\displaystyle(\mathcal{V}Tf)(z)=$ $\displaystyle\mathcal{V}D(zg(z))$
$\displaystyle=$ $\displaystyle zg(z)$ $\displaystyle=$ $\displaystyle
z(\mathcal{V}f)(z).$

So, $\mathcal{V}T=M_{z}\mathcal{V}$ and $\mathcal{V}T\mathcal{V}^{-1}=M_{z}$.
That is to say $\mathcal{V}$ transforms the operator $T$ into the operator
multiplication by $z$ on $\mathcal{S}_{0}^{2}(\mathbb{D})$. ∎

We can summarize the theorem by the following commutative diagram

$\mathcal{S}_{0}^{2}(\mathbb{D})$$\mathcal{S}_{0}^{2}(\mathbb{D})$$\mathcal{H}^{2}(\mathbb{D})$$\mathcal{H}^{2}(\mathbb{D})$$M_{z}$
$\mathcal{V}$ $\mathcal{V}$ $T$ $\mathcal{V}T=M_{z}\mathcal{V}$

######  Definition 2.

An element $a$ in Banach algebra $\mathcal{A}$ is called cyclic if the
subalgebra generated by $a$ is dense in $\mathcal{A}$.

The next proposition provides the relation between closed invariant subspaces
of multiplication operator and closed ideals of a Banach algebra. Proof of
this proposition can be found in [6].

######  Proposition 2.

Let $\mathcal{A}$ be a Banach algebra. Then the invariant subspaces of
multiplication by a cyclic element are exactly the closed ideals of
$\mathcal{A}$.

######  Lemma 1.

$J$ is a closed ideal of $\mathcal{S}_{0}^{2}(\mathbb{D})$ if and only if $J$
is an ideal of $\mathcal{S}^{2}(\mathbb{D})$ contained in
$\mathcal{S}_{0}^{2}(\mathbb{D})$.

###### Proof.

Using of Corollary 1, we see that for any $h\in\mathcal{S}^{2}(\mathbb{D})$
there exists $h_{1}\in\mathcal{S}_{0}^{2}(\mathbb{D})$ such that $h=c+h_{1}$.

Suppose $J$ is a closed ideal of $\mathcal{S}_{0}^{2}(\mathbb{D})$, then for
any $h\in\mathcal{S}^{2}(\mathbb{D})$ and $j\in J$,

$hj=(c+h_{1})j=cj+h_{1}j\in J.$

Since norm on both spaces are the same, $J$ is a closed ideal of
$\mathcal{S}^{2}(\mathbb{D})$. On the other hand, if $J$ is an ideal of
$\mathcal{S}^{2}(\mathbb{D})$ contained in $\mathcal{S}_{0}^{2}(\mathbb{D})$,
then it is clear that $J$ is a closed ideal of
$\mathcal{S}_{0}^{2}(\mathbb{D})$. ∎

######  Definition 3.

Let $K$ be a closed subset of the unit circle $\partial\mathbb{D}$. For and
inner function $G$ we say $G$ is associated with $K$ if

  * •

if $a_{1},a_{2},...$ are the zeros of $G(z)$ in the open disk, then all the
limit points of $\\{a_{k}\\}$ belong to $K$;

  * •

the measure determining the singular part of $G$ is supported on $K$.

The following theorem from Korenbljum [5] characterizes all closed ideals of
$\mathcal{S}^{2}(\mathbb{D})$. This theorem is essential in proving our main
theorem.

######  Theorem 2.

Suppose $K$ is a is closed subset of $\partial\mathbb{D}$ and let $G$ be an
inner function associated with $K$. Let $I_{\mathcal{S}^{2}(\mathbb{D})}(G;K)$
be the set of all $f\in\mathcal{S}^{2}(\mathbb{D})$ which are divisible by $G$
and which vanish on $K$. Then $I_{\mathcal{S}^{2}(\mathbb{D})}(G;K)$ is a
closed ideal of $\mathcal{S}^{2}(\mathbb{D})$. Moreover, every closed ideal of
$\mathcal{S}^{2}(\mathbb{D})$ is obtained in this manner.

######  Corollary 2.

Suppose $K$ is a is closed subset of $\partial\mathbb{D}$ and let $G$ be an
inner function associated with $K$. Let
$I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)$ be the set of all
$f\in\mathcal{S}_{0}^{2}(\mathbb{D})$ which are divisible by $G$ and which
vanish on $K$. Then $I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)$ is a closed
ideal of $\mathcal{S}_{0}^{2}(\mathbb{D})$. Moreover, every closed ideal of
$\mathcal{S}_{0}^{2}(\mathbb{D})$ is obtained in this manner.

###### Proof.

Clearly $I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)$ is a closed ideal of
$\mathcal{S}_{0}^{2}(\mathbb{D})$. Conversely, suppose $J$ be any ideal of
$\mathcal{S}_{0}^{2}(\mathbb{D})$. From the Lemma 1, $J$ is also an ideal of
$\mathcal{S}^{2}(\mathbb{D})$ contained in $\mathcal{S}_{0}^{2}(\mathbb{D})$.
By the Theorem 2, $J$ is of the form $I_{\mathcal{S}^{2}(\mathbb{D})}(G;K)$
for some $G$ and $K$ defined in Definition 3. Hence
$J=I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)$. ∎

######  Theorem 3 (Main Theorem).

Let $T$ be an operator

$(Tf)(z)=zf(z)+\int_{0}^{z}f(w)dw.$

defined on $\mathcal{H}^{2}(\mathbb{D})$. Then the lattice of closed invariant
subspaces is

$Lat~{}T=\\{S\subset\mathcal{H}^{2}(\mathbb{D}):S=\\{Df:f\in
I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)\\}~{}\text{for}~{}G,K~{}\text{defined~{}in~{}Definition
\ref{def3}}\\}.$

###### Proof.

From Corollary 1,

$\mathcal{S}_{0}^{2}(\mathbb{D})=\overline{span}\\{z^{n}:n\in\mathbb{N}\\},$

so $z$ is a cyclic element of the Banach algebra
$\mathcal{S}_{0}^{2}(\mathbb{D})$. Thus from the Proposition 2, closed
invariant subspaces of $M_{z}$ on $\mathcal{S}_{0}^{2}(\mathbb{D})$ are
exactly the closed ideals of $\mathcal{S}_{0}^{2}(\mathbb{D})$. Using
Corollary 2, the lattice of closed invariant subspace of $M_{z}$ acting on
$\mathcal{S}_{0}^{2}(\mathbb{D})$ is given by

$Lat~{}M_{z}=\\{I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K):G,K~{}\text{defined~{}in~{}Definition
\ref{def3}}\\}.$

Since $\mathcal{V}T\mathcal{V}^{-1}=M_{z}$, we see that
$\mathcal{V}^{-1}(I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K))$ is a closed
invariant subspace of $T$. From Theorem 1, we know that
$\mathcal{V}^{-1}(f)=Df$. So, $S=\\{Df:f\in
I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)\\}$ is a closed invariant subspace of
$T$. Hence,

$Lat~{}T=\\{S\subset\mathcal{H}^{2}(\mathbb{D}):S=\\{Df:f\in
I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)\\}~{}\text{for}~{}G,K~{}\text{defined~{}in~{}Definition
\ref{def3}}\\}.$

∎

######  Remark 2.

One can notice that if $K$ is of Lebesgue arc length measure zero, then the
ideal $I_{\mathcal{S}_{0}^{2}(\mathbb{D})}(G;K)\\}$ is the zero ideal.

## References

  * [1] A. Aleman and B. Korenblum, _Volterra invariant subspaces of $H^{p}$,_ Bull. Sci. Math. 132 (2008).
  * [2] A. Beurling, _On two problems concerning linear transformations in Hilbert space,_ Acta Math. 81, (1948).
  * [3] C. C. Cowen, G. Gunatillake and E. Ko, _Hermitian weighted composition operators and Bergman extremal functions_ , Complex Anal. Oper. Theory 7 (2013), no. 1, 69–99.
  * [4] Y. Katznelson, _An introduction to harmonic analysis_ , Third edition, Cambridge Mathematical Library, Cambridge University Press, Cambridge, (2004).
  * [5] B. I. Korenbljum, _Invariant subspaces of the shift operator in a weighted Hilbert space,_ (Russian) Mat. Sb. (N.S.) 89(131) (1972).
  * [6] A. Montes-Rodriguez, M. Ponce-Escudero and S. A. Shkarin, _Invariant subspaces of parabolic self-maps in the Hardy space,_ Math. Res. Lett. 17 (2010).
  * [7] W. Rudin, _The closed ideals in an algebra of analytic functions_ , Canad. J. Math. 9 (1957), 426-434.
  * [8] D. Sarason, _A remark on the Volterra operator,_ J. Math. Anal. Appl. 12 (1965).
  * [9] D. Sarason, _Invariant subspaces,_ Topics in operator theory, Math. Surveys, No. 13, Amer. Math. Soc., Providence, R.I., (1974).

Department of Mathematics and Statistics, The University of Toledo, Toledo,
Ohio, E-mail address<EMAIL_ADDRESS>
Department of Mathematics and Statistics, The University of Toledo, Toledo,
Ohio, Email address<EMAIL_ADDRESS>


# On the moments of the Ulam-Kac adder

Gage Bonner<EMAIL_ADDRESS>Department of Mathematics, University of
Wisconsin – Madison, WI 53706, USA

###### Abstract.

Let $\\{U(n)\\}_{n\geq 0}$ be a sequence of independent random variables such
that $U(n)$ is distributed uniformly on $\\{0,1,2\dots n\\}$. The Ulam-Kac
adder is the history-dependent random sequence defined by
$X_{n+1}=X_{n}+X_{U(n)}$ with the initial condition $X_{0}=1$. We show that
for each $m\geq 1$, it holds that $\log E[X_{n}^{m}]/\sqrt{n}$ approaches a
constant $c_{m}$ as $n\to\infty$. Loose bounds are provided for the constants
$c_{m}$.

##  1\. Introduction

History-dependent random sequences are sequences in which the distribution of
the outcome of the $n^{\text{th}}$ step depends non-trivially on the entire
history of the process. Such sequences have numerous applications in sciences
such as in polymer physics where basic models of confined polymer chains in
small enclosures consist of lattice sites which are occupied sequentially. The
points evolve according to a self-avoiding random walk [8, 11, 17] and the
investigation of the analytical properties of such processes remains an active
area of research [5, 18].

The focus of this paper is the so-called Ulam-Kac adder, following the “Ulam-
Kac process” defined in [7]. Let $\\{U(n)\\}_{n\geq 0}$ be a sequence of
independent random variables such that $U(n)$ is distributed uniformly on
$\\{0,1,2\dots n\\}$. The Ulam-Kac adder is the history-dependent random
sequence defined by

$X_{n+1}=X_{n}+X_{U(n)};\quad X_{0}=1.$ (1.1)

This sequence was first investigated by Ulam [26]. The path space for the
first five steps of this sequence is shown in Figure 1.

$1$$2$$3$$4$$4$$5$$6$$8$$6$$5$$5$$6$$7$$8$$6$$7$$8$$10$$7$$8$$9$$12$$16$$12$$10$$9$$12$$10$$8$$7$$10$$9$$7$$6$
Figure 1. Illustration of the paths through the Ulam-Kac adder for the first
five steps. Each of the paths are realized with equal probability.

The first explicit results on the behavior of the sequence were obtained by
Kac [14] where it was shown using generating functions and the method of
steepest descent that

$\displaystyle\log E[X_{n}]$ $\displaystyle\sim 2\sqrt{n},$ (1.2)
$\displaystyle\log E[X_{n}^{2}]$
$\displaystyle\sim\sqrt{2n\mathopen{}\mathclose{{}\left(5+\sqrt{17}}\right)}.$
(1.3)

where we write $f(n)\sim g(n)$ if $\lim_{n\to\infty}f(n)/g(n)=1$. In [16], Eq.
(1.1) is generalized to

$X_{n+1}=X_{n}+\gamma X_{U(n)};\quad X_{0}=1,$ (1.4)

and the analogues of Eqs. (1.2) and (1.3) are obtained. The Ulam-Kac adder was
also considered by [3] in the context of extensions of random Fibonacci
sequences. There, the following results are conjectured based on numerical
evidence:

######  Conjecture 1.1.

There is an increasing sequence of positive constants $(c_{m})_{m\geq 1}$ such
that

$\log E[X_{n}^{m}]\sim c_{m}\sqrt{n}\quad\forall m\geq 1.$ (1.5)

######  Conjecture 1.2.

There are positive constants $\mu,\sigma$ such that

$\frac{\log X_{n}-\mu n^{1/2}}{\sigma
n^{1/4}}\to\mathcal{N}(\mu,\sigma^{2})\text{ in law}.$ (1.6)

In [7], a Poisson-regularized continuous time analogue of the Ulam-Kac adder
was introduced as an alternate method of analyzing the asymptotic behavior of
moments. This was applied to the first and second moment in several
generalizations and also in [21] to analyze the tail distribution of related
sequences.

We mention here a connection between the Ulam-Kac adder and addition chains
which does not appear to have been explicitly stated previously. An addition
chain for $n$ is a sequence $(a_{0},a_{1},a_{2}\dots a_{m-1},a_{m})$ such that
$a_{0}=1,a_{m}=n$ where $a_{i}=a_{j}+a_{k}$ for each $1\leq i\leq m$ and some
$0\leq j,k<i$. An open problem in computer science is to obtain efficiently
the shortest addition chain for arbitrary $n$ [22, 6, 24, 23]. This remains an
active area of research in high-performance computing [2, 15]. In the case
$j=i-1$, the chain is said to have taken a _star step._ With reference to
Figure 1, we note that each realization of the Ulam-Kac adder is addition
chain consisting only of star steps. Such chains are called Brauer chains, or
simply star chains. Computing $P(X_{k}=n)$ is equivalent to computing the
number of Brauer chains for $n$ that have length $k$. The first passage times
of the Ulam-Kac adder are therefore directly related to minimal length Brauer
chains.

The main result of this paper is Theorem 1.3, which provides an affirmative
answer to Conjecture 1.1 as well as bounds on the constants $c_{m}$ appearing
therein.

######  Theorem 1.3 (Main result).

For all $m\geq 1$, there exists a
$c_{m}\in[2m,2^{\mathopen{}\mathclose{{}\left(\tfrac{m}{2}+1}\right)e^{3\sqrt{m}}}]$
such that

$\log E[X_{n}^{m}]\sim c_{m}\sqrt{n}.$ (1.7)

The remainder of this paper is organized as follows. In Section 2, we provide
a brief review the basic method of computation for the $m=1,2$ moments. Then,
Section 3 contains our outline of the main steps of the proof. In Section 4 we
apply our main result to compute $E[X_{n}^{3}]$ explicitly and compare the its
value to a previous numerical estimation. Section 5 contains all the proofs of
the intermediate results required to obtain the main theorem. Finally, Section
6 contains our remarks on the tightness of the bounds in the main result as
well as some additional conjectures related to this work. Appendix A contains
supplementary proofs.

##  2\. The $m=1,2$ moments

The basic method for the analysis of the Ulam-Kac adder is to condition on
possible values of $U(n)$ and telescoping the resulting series [14], which we
will now demonstrate. Defining $\mu_{n}\coloneqq E[X_{n}]$, Eq. (1.1) gives

$\mu_{n+1}=\mu_{n}+E[X_{U(n)}]=\mu_{n}+\frac{1}{n+1}\sum_{\ell=0}^{n}\mu_{\ell},$
(2.1)

where we have applied $E[X_{U(n)}|U(n)=\ell]=\mu_{\ell}.$ We then take $n\to
n+1$ in Eq. (2.1) and eliminate the remaining sum to obtain

$\displaystyle\mu_{n+2}$
$\displaystyle=\mu_{n+1}+\frac{1}{n+2}\mathopen{}\mathclose{{}\left[mu_{n+1}+\sum_{m=0}^{n}\mu_{m}}\right]$
(2.2) $\displaystyle 0$
$\displaystyle=(n+2)\mu_{n+2}-2(n+2)\mu_{n+1}+(n+1)\mu_{n},$ (2.3)

along with the initial conditions $\mu_{0}=1,\mu_{1}=2$. The classical
Laguerre polynomials $L_{n}(x)$ satisfy the recurrence relation [25]

$(2+k)L_{k+2}(x)-(2k+3-x)L_{k+1}(x)+(k+1)L_{k}(x)=0,\quad L_{0}=1,L_{1}=1-x.$
(2.4)

Comparing this to Eq. (2.3), we see that they coincide when $x=-1$ and
therefore

$\mu_{n}=L_{n}(-1)=\sum_{\ell=0}^{n}\binom{n}{\ell}\frac{1}{\ell!}.$ (2.5)

This directly leads to $\log\mu_{n}\sim 2\sqrt{n}$ using the known asymptotics
of Laguerre polynomials. See [3, 16] for alternative calculations of the same
result.

We will now compute the second moment of the Ulam-Kac adder
$\sigma_{n}\coloneqq E[X_{n}^{2}]$ by determining an ODE satisfied by the
generating function of the sequence $\\{\sigma_{n}\\}_{n\geq 0}$. By squaring
Eq. (1.1), we obtain

$X_{n+1}^{2}=X_{n}^{2}+2X_{n}X_{U(n)}+X_{U(n)}^{2}.$ (2.6)

We would like to calculate the expected value of this equation, but the cross
term $X_{n}X_{U(n)}$ is not easily identified with $\sigma_{n}$. To obtain a
closed system of recurrence relations, we define
$\alpha_{n}\coloneqq\sum_{\ell=0}^{n-1}E[X_{n}X_{\ell}]$, then taking the
expectation of Eq. (2.6) gives

$\sigma_{n+1}=\sigma_{n}+\frac{2}{n+1}\mathopen{}\mathclose{{}\left(\sigma_{n}+\alpha_{n}}\right)+\frac{1}{n+1}\sum_{\ell=0}^{n}\sigma_{\ell}.$
(2.7)

Furthermore, by applying Eq. (1.1) to $\alpha_{n+1}$, we obtain

$\alpha_{n+1}=\sigma_{n}+\alpha_{n}+\frac{1}{n+1}\sum_{\ell=0}^{n}\mathopen{}\mathclose{{}\left(\sigma_{\ell}+2\alpha_{\ell}}\right).$
(2.8)

We define the generating functions $G(z)=\sum_{n=0}^{\infty}\sigma_{n}z^{n}$
and $M(z)=\sum_{n=0}^{\infty}\alpha_{n}z^{n}$ which have radii of convergence
equal to 1 as will be shown in Section 3. Multiplying Eqs. (2.7) and (2.8)
through by $(n+1)z^{n}$ and summing the resulting expressions over $n\geq 0$
gives

$\displaystyle G^{\prime}$
$\displaystyle=(zG)^{\prime}+2(G+M)+\frac{1}{1-z}G,\quad G(0)=1$ (2.9a)
$\displaystyle M^{\prime}$
$\displaystyle=(zG)^{\prime}+(zM)^{\prime}+\frac{1}{1-z}(G+2M),\quad M(0)=0.$
(2.9b)

We can eliminate $M(z)$ using Eq. (2.9a) to obtain

$(z^{2}-6z+7)G(z)+(3z-8)(1-z)^{2}G^{\prime}(z)+(1-z)^{4}G^{\prime\prime}(z)=0,\quad
G(0)=1,G^{\prime}(0)=4,$ (2.10)

which is the ODE we sought. Substitution of the asymptotic form $e^{S(x)}$
into Eq. (2.10) and application of the method of dominant balance [4] provides
that
$G(z)\sim\exp\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(\frac{5+\sqrt{17}}{2}}\right)\frac{1}{1-z}}\right]$
near $z=1$ which implies that $\log\sigma_{n}\sim\sqrt{2(5+\sqrt{17})n}$.

##  3\. Outline of Proof

We will begin by constructing a closed set of recurrence relations for the
$m^{\text{th}}$ moment. Then, we will create the system of linear differential
equations for the associated generating functions. By closely studying the
properties of the matrices defining this system, we are able to apply certain
technical results for the asymptotic behavior of their solutions and thereby
obtain Theorem 1.3.

Let $\mathcal{M}_{\geq 1}$ be the set which contains the empty set $\emptyset$
plus all nonempty multisets of finite cardinality on positive integers. We
write an element $\mathfrak{a}\in\mathcal{M}_{\geq 1}$ as
$\mathfrak{a}=[p_{1},p_{2},p_{3}\dots p_{n}]$, where $p_{i}\in\mathbb{Z}_{\geq
1}$ for all $1\leq i\leq n$ and it is allowed that $p_{i}=p_{j}$ for $i\neq
j$. We call each $p_{i}$ a _part_ of $\mathfrak{a}$ and write
$|\mathfrak{a}|=n$ when $\mathfrak{a}$ has $n$ parts. We reserve
$\lVert\mathfrak{a}\rVert$ to denote the sum of all parts of $|\mathfrak{a}|$.
We also recall that addition and subtraction can be defined on multisets as
follows. If $|\mathfrak{a}|,|\mathfrak{b}|\in\mathcal{M}_{\geq 1}$ where
$|\mathfrak{a}|=m$ and $|\mathfrak{b}|=n$ then $\mathfrak{a}+\mathfrak{b}$ is
the multiset with $m+n$ parts created from the disjoint union of all parts of
$\mathfrak{a}$ and $\mathfrak{b}$. Also, $\mathfrak{a}-\mathfrak{b}$ is the
unique multisubset $\mathfrak{c}\subseteq\mathfrak{a}$ such that
$\mathfrak{c}+\mathfrak{b}=\mathfrak{a}$. Successive addition and subtraction
should be applied from left to right, e.g. we have
$\mathfrak{a}+\mathfrak{b}-\mathfrak{c}=(\mathfrak{a}+\mathfrak{b})-\mathfrak{c}$.

Let $\mathcal{P}(n)=\\{\mathfrak{p}\in\mathcal{M}_{\geq
1}:\lVert\mathfrak{p}\rVert=n\\}$ be the set of all integer partitions of
$n\geq 0$ with the definition $\mathcal{P}(0)=\emptyset$. We define
$\mathcal{C}^{(m)}$ as the set of all terms that can occur in our system of
recurrence relations for the $m^{\text{th}}$ moment,

$\mathcal{C}^{(m)}=\\{C_{n}^{(m)}(q;\mathfrak{p}):n\geq 1,1\leq q\leq
m,\mathfrak{p}\in\mathcal{P}(m-q)\\},\quad m\geq 1,$ (3.1)

such that each element of $\mathcal{C}^{(m)}$ is given by

$C_{n}^{(m)}(q;\mathfrak{p})=\sum_{\alpha_{1}=0}^{n-1}\sum_{\alpha_{2}=0}^{n-1}\cdots\sum_{\alpha_{|\mathfrak{p}|}=0}^{n-1}E\mathopen{}\mathclose{{}\left[X_{n}^{q}\prod_{\ell=1}^{|\mathfrak{p}|}X^{p_{\ell}}_{\alpha_{\ell}}}\right].$
(3.2)

Motivated by the patterns observed in Eqs. (2.7) and (2.8), we generally
expect to obtain recurrence relations each with three types of terms:
$x_{n},x_{n}/(n+1)$ and $\tfrac{1}{n+1}\sum_{\ell=0}^{n}x_{\ell}$. We
therefore introduce the following notation:

$\displaystyle\overline{x}_{n}$ $\displaystyle\coloneqq\frac{x_{n}}{n+1},$
(3.3a) $\displaystyle\hat{x}_{n}$
$\displaystyle\coloneqq\frac{1}{n+1}\sum_{\ell=0}^{n}x_{\ell}.$ (3.3b)

We can now compute the most general recurrence relation, which is the result
of the following Lemma whose proof is found in Section 5.1.

######  Lemma 3.1.

For a fixed $m\geq 1$, each
$C_{n+1}^{(m)}(q;\mathfrak{p})\in\mathcal{C}^{(m)}$ satisfies

$\displaystyle C_{n+1}^{(m)}(q;\mathfrak{p})$
$\displaystyle=\sum_{\mathfrak{b}\subseteq\mathfrak{p}}C_{n}^{(m)}\mathopen{}\mathclose{{}\left(q+\lVert\mathfrak{b}\rVert;\mathfrak{p}-\mathfrak{b}}\right)$
(3.4a)
$\displaystyle+\sum_{\beta=1}^{q-1}\binom{q}{\beta}\sum_{\mathfrak{b}\subseteq\mathfrak{p}+[q-\beta]}\overline{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(\beta+\lVert\mathfrak{b}\rVert;\mathfrak{p}+[q-\beta]-\mathfrak{b}}\right)$
(3.4b)
$\displaystyle+\sum_{\begin{subarray}{c}\mathfrak{b}\subseteq\mathfrak{p}+[q]\\\
\mathfrak{b}\neq\emptyset\end{subarray}}\hat{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(\lVert\mathfrak{b}\rVert;\mathfrak{p}+[q]-\mathfrak{b}}\right).$
(3.4c)

We will show how Eq. (3.4) can be written as a finite linear system of first
order difference equations. This will allow us to obtain preliminary results
concerning the asymptotic behavior of the system. We first collect the
$C_{n}^{(m)}(q;\mathfrak{p})$ into a vector, $h^{(m)}_{n}$. The following
definition is important for this construction. For a fixed $m\geq 1$, the
vector $h^{(m)}_{n}=(C_{n}^{(m)}(q;\mathfrak{p}))$ is said to be in _canonical
ordering_ if its entries are sorted first by ascending $|\mathfrak{p}|$ then
by descending $q$. The number of entries in this vector is $\mathfrak{C}(m)$.

Now, Eq. (3.4) can first be expressed generically as

$h^{(m)}_{n+1}=\mathbf{C}_{1}h^{(m)}_{n}+\frac{1}{n+1}\mathbf{C}_{2}h^{(m)}_{n}+\frac{1}{n+1}\sum_{\ell=0}^{n}\mathbf{C}_{3}h^{(m)}_{\ell}.$
(3.5)

Each of the matrices $\mathbf{C}_{1},\mathbf{C}_{2}$ and $\mathbf{C}_{3}$ have
entries which depend only on $m$ and can be extracted from Eqs. (3.4a), (3.4b)
and (3.4c), respectively. Define

$u^{(m)}_{n+1}=\sum_{\ell=0}^{n}\mathbf{C}_{3}h^{(m)}_{\ell}=\mathbf{C}_{3}h^{(m)}_{n}+u^{(m)}_{n}.$
(3.6)

Then, we have

$\begin{pmatrix}h^{(m)}_{n+1}\\\
u^{(m)}_{n+1}\end{pmatrix}=\begin{pmatrix}\mathbf{C}_{1}+\tfrac{1}{n+1}\mathbf{C}_{2}&\tfrac{1}{n+1}\mathbf{C}_{3}\\\
\mathbf{C}_{3}&\mathbf{I}\end{pmatrix}\begin{pmatrix}h^{(m)}_{n}\\\
u^{(m)}_{n}\end{pmatrix}.$ (3.7)

We can now apply the following result due to Pituk on asymptotically constant
systems of linear difference equations.

######  Theorem 3.2 ([20], Theorem 1).

Let $x=(x_{n})_{n\in\mathbb{N}}$ be a solution of the difference equation
$x_{n+1}=\mathbf{L}_{n}x_{n},n\in\mathbb{N}$. If
$\lim_{n\to\infty}\mathbf{L}_{n}=\mathbf{L}$ for a constant matrix
$\mathbf{L}$ then either $x_{n}=0$ for all large $n$ or the limit

$\rho(x)=\lim_{n\to\infty}||x_{n}||^{1/n}$ (3.8)

exists and is equal to the modulus of one of the eigenvalues of $\mathbf{L}$.

Examining Eq. (3.7), Theorem 3.2 applies with

$\mathbf{L}=\begin{pmatrix}\mathbf{C}_{1}&\mathbf{0}\\\
\mathbf{C}_{3}&\mathbf{I}\end{pmatrix}.$ (3.9)

We require the following proposition to compute the eigenvalues of
$\mathbf{L}$.

######  Proposition 3.3.

The matrix $\mathbf{C_{1}}$ is lower triangular and each entry on its diagonal
is $1$.

###### Proof.

The entries of $\mathbf{C}_{1}$ are controlled by Eq. (3.4a); it implies that
$(q;\mathfrak{p})$ depends on each
$(q+\lVert\mathfrak{b}\rVert;\mathfrak{p}-\mathfrak{b})$ for
$\mathfrak{b}\subseteq\mathfrak{p}$. Each
$(q+\lVert\mathfrak{b}\rVert;\mathfrak{p}-\mathfrak{b})$ comes before
$(q;\mathfrak{p})$ in the canonical ordering, which implies the lower
triangular structure. There is a unique $\mathfrak{b}$ such that
$(q+\lVert\mathfrak{b}\rVert;\mathfrak{p}-\mathfrak{b})=(q;\mathfrak{p})$,
namely the empty set, and hence there are ones on the diagonal of
$\mathbf{C}_{1}$. ∎

By Proposition 3.3, $\mathbf{L}$ and by Eq. (3.9), $\mathbf{L}$ is a lower
triangular matrix whose diagonal entries are all equal to 1, and hence
$\mathbf{L}$ has one eigenvalue equal to 1 of multiplicity $2\mathfrak{C}(m)$.
Hence, by Theorem 3.2, we have that each increasing solution $x_{n}$ of Eq.
(3.7) satisfies

$\lim_{n\to\infty}||x_{n}||^{1/n}=1.$ (3.10)

We remark that a theorem of Birkhoff and Trjitzinsky (Theorem 1 of [27]) then
implies that $\log(x_{n})\sim n^{p}$ for a rational $p\in(0,1)$.

To each sequence $(C_{n}^{(m)}(q;\mathfrak{p}))_{n\geq 0}$ we associate a
generating function $G^{(m)}_{q,\mathfrak{p}}(z)$ defined by

$G^{(m)}_{q,\mathfrak{p}}(z)=\sum_{n=0}^{\infty}C_{n}^{(m)}(q;\mathfrak{p})z^{n},$
(3.11)

which has radius of convergence equal to 1 since the growth of
$(C_{n}^{(m)}(q;\mathfrak{p}))_{n\geq 0}$ is sub-exponential. We will apply
Lemma 3.1 to obtain a system of linear differential equations satisfied by
these generating functions. Therefore, we collect them into a vector
$g^{(m)}(z)=(G^{(m)}_{q,\mathfrak{p}}(z))$ in the canonical ordering. We will
refer to generating functions by just their indices $(q;\mathfrak{p})$. Also,
$(q;\mathfrak{p})_{i}$ will refer to the $i^{\text{th}}$ generating function
when listed in the canonical ordering. We will similarly refer to the
partition appearing in the label of the $i^{\text{th}}$ generating function as
$\mathfrak{p}_{i}$. The $m$ label will generally be omitted for brevity except
where required.

Now, we multiply each term in Eq. (3.4) by $(n+1)z^{n}$ and sum the resulting
equation over all $n\geq 0$ to transform the recurrence relations to
differential equations. Recalling that the generating functions of
$a_{n},\overline{a}_{n}$ and $\hat{a}_{n}$ are $(zA(z))^{\prime},A(z)$ and
$A(z)/(1-z)$, respectively, we obtain the following system of differential
equations

$g^{\prime}(z)=z\mathbf{C}_{1}g^{\prime}(z)+\mathopen{}\mathclose{{}\left[\mathbf{C}_{1}+\mathbf{C}_{2}+\frac{1}{1-z}\mathbf{C}_{3}}\right]g(z).$
(3.12)

The matrix $(\mathbf{I}-z\mathbf{C}_{1})^{-1}$ is invertible for $0\leq z<1$
(by Proposition 3.3), so we can write

$g^{\prime}(z)=\mathbf{R}(z)g(z),\quad\mathbf{R}(z)=(\mathbf{I}-z\mathbf{C}_{1})^{-1}\mathopen{}\mathclose{{}\left[\mathbf{C}_{1}+\mathbf{C}_{2}+\frac{1}{1-z}\mathbf{C}_{3}}\right].$
(3.13)

We make the substitution

$g(z)=\mathbf{P}(z)h(z),\quad\mathbf{P}(z)=\text{diag}\\{(1-z)^{-|\mathfrak{p}_{1}|},(1-z)^{-|\mathfrak{p}_{2}|},\dots,(1-z)^{-|\mathfrak{p}_{\mathfrak{C}(m)}|}\\},$
(3.14)

as well as $u=1/(1-z)$ which results in a new system

$h^{\prime}(u)=\mathbf{S}(u)h(u).$ (3.15)

In Section 5.2, we establish the following lemma which collects the relevant
properties of $\mathbf{S}(u)$.

######  Lemma 3.4.

For $u\geq 1$, and for some $N\geq 1$ there is a sequence of matrices
$\\{\mathbf{M}_{i}\\}_{0\leq i\leq N}$ such that $\mathbf{M}_{0}$ is primitive
with Perron-Frobenius eigenvalue bounded above by
$2^{\mathopen{}\mathclose{{}\left(\tfrac{m}{2}+1}\right)e^{3\sqrt{m}}}$ and
such that the matrix $\mathbf{S}(u)$ defined in Eq. (3.15) can be expressed as

$\mathbf{S}(u)=\sum_{i=0}^{N}u^{-i}\mathbf{M}_{i}.$ (3.16)

To obtain the asymptotic behavior of $h(u)$, we will apply the following
lemma, a consequence of the Hartman-Wintner theorem [13] whose proof is in
Section 5.3.

######  Lemma 3.5.

Let $d,n\geq 1$, let $(\mathbf{M}_{i})_{0\leq i\leq n}$ be $d\times d$
matrices for each $i$ and let $y(t)=(y_{i})_{1\leq i\leq d}$ be a column
vector. Consider the system of differential equations

$y^{\prime}(t)=\mathopen{}\mathclose{{}\left(\sum_{\ell=0}^{n}\mathbf{M}_{\ell}t^{-\ell}}\right)y(t);\quad
y(1)=y_{0},t\geq 1.$ (3.17)

Then, if $\mathbf{M}_{0}$ is a non-negative, irreducible matrix with Perron-
Frobenius eigenvalue $\lambda$, it holds for some $\alpha,\beta>0$ that

$y_{i}(t)\sim\alpha(t^{\beta}+o(t^{\beta}))e^{\lambda t},\quad 1\leq i\leq d.$
(3.18)

We can now prove Theorem 1.3.

###### Proof of Theorem 1.3.

By Lemmas 3.4 and 3.5, for a given $m\geq 1$, Eq. (3.15) admits solutions
$(h_{i}^{(m)}(u))_{1\leq i\leq\mathfrak{C}(m)}$ such that for some
$\alpha,\beta>0$,

$h_{i}^{(m)}(u)\sim\alpha(u^{\beta}+o(u^{\beta}))e^{\lambda_{m}u},\quad 1\leq
i\leq\mathfrak{C}(m),$ (3.19)

where $\lambda_{m}$ is the Perron-Frobenius eigenvalue of the matrix
$\mathbf{M}_{0}$ defined in Lemma 3.4. Hence, we have
$G_{m,\emptyset}^{(m)}(z)\sim\frac{\alpha}{(1-z)^{\beta}}\exp\mathopen{}\mathclose{{}\left(\frac{\lambda_{m}}{1-z}}\right)$
which by the method of saddle point asymptotics [10] implies that $\log
E[X_{n}^{m}]\sim 2\sqrt{\lambda_{m}n}$.

We will now put the required bounds on $\lambda_{m}$. Note that a lower bound
can be obtained quickly; by Jensen’s inequality, we have

$E[X_{n}^{m}]\geq E[X_{n}]^{m}\implies\log E[X_{n}^{m}]\sim
c_{m}\sqrt{n}\quad\text{for }c_{m}\geq 2m.$ (3.20)

To obtain an upper bound, we apply Lemma 3.4 to obtain $\lambda_{m}\leq
2^{\mathopen{}\mathclose{{}\left(\tfrac{m}{2}+1}\right)e^{3\sqrt{m}}}$ and
hence that

$\log E[X_{n}^{m}]\sim c_{m}\sqrt{n}\quad\text{for }c_{m}\leq
2^{\mathopen{}\mathclose{{}\left(\tfrac{m}{2}+1}\right)e^{3\sqrt{m}}}.$ (3.21)

Combining Eqs. (3.20) and (3.21) completes the proof. ∎

##  4\. Example: $m=3$

We will now compute the exact value of the third moment of the Ulam-Kac adder,
$E[X_{n}^{3}]$. This has not been explicitly computed previously, although [3]
estimated $\log E[X_{n}^{3}]\sim 6.5\sqrt{n}$ based on numerical evidence.
When $m=3$, there are four types of sequences in $\mathcal{C}^{(3)}$, namely

$\displaystyle\gamma_{n}$ $\displaystyle\coloneqq
C_{n}^{(3)}(3,\emptyset)=E[X_{n}^{3}],$ (4.1a) $\displaystyle a_{n}$
$\displaystyle\coloneqq
C_{n}^{(3)}(2,\\{1\\})=\sum_{\ell=0}^{n-1}E[X_{n}^{2}X_{\ell}],$ (4.1b)
$\displaystyle b_{n}$ $\displaystyle\coloneqq
C_{n}^{(3)}(1,\\{2\\})=\sum_{\ell=0}^{n-1}E[X_{n}X_{\ell}^{2}],$ (4.1c)
$\displaystyle c_{n}$ $\displaystyle\coloneqq
C_{n}^{(3)}(1,\\{1,1\\})=\sum_{\ell=0}^{n-1}\sum_{m=0}^{n-1}E[X_{n}X_{\ell}X_{m}].$
(4.1d)

Applying Lemma 3.1, we obtain the following system of recurrence relations.

$\displaystyle\gamma_{n+1}$
$\displaystyle=\gamma_{n}+3\mathopen{}\mathclose{{}\left(2\overline{\gamma}_{n}+\overline{a}_{n}+\overline{b}_{n}}\right)+\hat{\gamma}_{n},$
(4.2a) $\displaystyle a_{n+1}$
$\displaystyle=\gamma_{n}+a_{n}+2\mathopen{}\mathclose{{}\left(\overline{\gamma}_{n}+2\overline{a}_{n}+\overline{c}_{n}}\right)+\mathopen{}\mathclose{{}\left(\hat{\gamma}_{n}+\hat{a}_{n}+\hat{b}_{n}}\right),$
(4.2b) $\displaystyle b_{n+1}$
$\displaystyle=\gamma_{n}+b_{n}+\mathopen{}\mathclose{{}\left(\hat{\gamma}_{n}+\hat{a}_{n}+\hat{b}_{n}}\right),$
(4.2c) $\displaystyle c_{n+1}$
$\displaystyle=\gamma_{n}+c_{n}+2a_{n}+\mathopen{}\mathclose{{}\left(\hat{\gamma}_{n}+3\hat{a}_{n}+3\hat{c}_{n}}\right).$
(4.2d)

We elaborate on how Eq. (4.2b) is obtained; the other equations follow
similarly. To calculate $c_{n}$ we set $m=3,q=2$ and $\mathfrak{p}=\\{1\\}$ in
Eq. (3.4). The set of subsets of $\mathfrak{p}$ is $\\{\\{1\\},\emptyset\\}$,
hence the term on the righthand side of Eq. (3.4a) is

$C_{n}^{(m)}(2+1;\emptyset)+C_{n}^{(m)}(2+0;\\{1\\})=\gamma_{n}+a_{n}.$ (4.3)

Consider now the term on Eq. (3.4b). Since $q=2$, the sum over $\beta$
contains only the term $\beta=1$, hence there will be an overall factor of
$\binom{2}{1}=2$. The sum over $\mathfrak{b}$ is a sum over the four subsets
of $\mathfrak{p}+[1]=\\{1,1\\}$, hence the entire term is

$2[\overline{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(1+2;\emptyset}\right)+2\overline{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(1+1;\\{1\\}}\right)+\overline{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(1+0;\\{1,1\\}}\right)]=2(\overline{\gamma}_{n}+2\overline{a}_{n}+\overline{c}_{n}).$
(4.4)

Finally, the term on Eq. (3.4c) is a sum over the three non-empty subsets of
$\mathfrak{p}+[2]=\\{1,2\\}$, giving

$\hat{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(3;\emptyset}\right)+\hat{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(2;\\{1\\}}\right)+\hat{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(1;\\{2\\}}\right)=\hat{\gamma}_{n}+\hat{a}_{n}+\hat{b}_{n}.$
(4.5)

Collecting all these terms gives Eq. (4.2b) as required. Next, we construct
the vector of generating functions in canonical order,
$g=(G^{(3)}_{3,\emptyset},G^{(3)}_{2,\\{1\\}},G^{(3)}_{1,\\{2\\}},G^{(3)}_{1,\\{1,1\\}})$.
Multiplying Eqs. (4.2) by $(n+1)z^{n}$ and summing over $n$ results in a
system of differential equations $g^{\prime}(z)=\mathbf{R}(z)g(z)$ where

$\mathbf{R}(z)=\frac{1}{(1-z)^{2}}\mathopen{}\mathclose{{}\left(\begin{array}[]{cccc}8-7z&3-3z&3-3z&0\\\
4z+\frac{1}{1-z}+3&6-2z&3z+1&2-2z\\\ 6z+\frac{1}{1-z}+1&3z+1&2(z+1)&0\\\
\frac{2-(z-1)z(2z+11)}{(z-1)^{2}}&5z-\frac{8}{z-1}-3&-\frac{z(3z+5)}{z-1}&3z+4\\\
\end{array}}\right).$ (4.6)

We make the transformation $g(z)=\mathbf{P}(z)h(z)$ where

$\mathbf{P}=\text{diag}\\{(1-z)^{-0},(1-z)^{-1},(1-z)^{-1},(1-z)^{-2}\\},$
(4.7)

as well as the substitution $u=1/(1-z)$, which results in the system

$h(u)^{\prime}=\mathopen{}\mathclose{{}\left(\mathbf{M}_{0}+\mathbf{M}_{1}u^{-1}+\mathbf{M}_{2}u^{-2}+\mathbf{M}_{3}u^{-3}}\right)h(u),$
(4.8)

where

$\mathbf{M}_{0}=\mathopen{}\mathclose{{}\left(\begin{array}[]{cccc}1&3&3&0\\\
1&4&4&2\\\ 1&4&4&0\\\ 2&8&8&7\\\ \end{array}}\right).$ (4.9)

Observe that $\mathbf{M}_{0}$ is a primitive matrix with Perron-Frobenius
eigenvalue $\lambda^{(3)}$ given by

$\displaystyle\lambda^{(3)}$
$\displaystyle=\max_{x\in\mathbb{R}}\\{x^{4}-16x^{3}+49x^{2}-10x=0\\}$ (4.10a)
$\displaystyle=\frac{2}{3}\mathopen{}\mathclose{{}\left(8+\sqrt{109}\cos\mathopen{}\mathclose{{}\left(\frac{1}{3}\tan^{-1}\mathopen{}\mathclose{{}\left(\frac{6\sqrt{22245}}{703}}\right)}\right)}\right)$
(4.10b) $\displaystyle=11.979293127704\dots$ (4.10c)

By the same arguments in the proof of Theorem 1.3, we conclude that

$\log E[X_{n}^{3}]\sim 2\sqrt{\lambda^{(3)}n}\approx 6.92\sqrt{n}.$ (4.11)

To compare this with the previous numerical estimate $\log E[X_{n}^{3}]\sim
6.5\sqrt{n}$ from [3], we reduce $g^{\prime}(z)=\mathbf{R}(z)g(z)$ to a single
differential equation for $\Gamma(z)\coloneqq G_{3,\emptyset}^{(3)}(z)$. We
find the following fourth-order differential equation,

$\sum_{n=0}^{4}f_{n}(z)\frac{\text{d}^{n}\Gamma}{\text{d}z^{n}}=0,\quad\Gamma(0)=1,\Gamma^{\prime}(0)=8,\Gamma^{\prime\prime}(0)=91,\Gamma^{\prime\prime\prime}(0)=1258,$
(4.12)

where

$\displaystyle f_{0}(z)$ $\displaystyle=17-44z+21z^{2}-2z^{3},$ $\displaystyle
f_{1}(z)$ $\displaystyle=(1-z)(-236+377z-173z^{2}+22z^{3}),$ $\displaystyle
f_{2}(z)$ $\displaystyle=(1-z)^{3}(174-156z+31z^{2}),$ $\displaystyle
f_{3}(z)$ $\displaystyle=(1-z)^{5}(-27+11z),$ $\displaystyle f_{4}(z)$
$\displaystyle=(1-z)^{7}.$

We can construct a power series solution with $N$ terms to Eq. (4.12) to
calculate the exact values of the sequence $(E[X_{n}^{3}])_{0\leq n\leq N}$
and thus obtain numerical estimates for its asymptotic behavior. Substituting
$\Gamma(z)=\sum_{n=0}^{\infty}\gamma_{n}z^{n}$ into Eq. (4.12) gives the
following recurrence relation for $(\gamma_{n})_{n\geq 0}$,

$\displaystyle\sum_{i=0}^{7}g_{i}(n)\gamma_{n+i}=0,\quad\text{where}$
$\displaystyle\gamma_{0}=1,\gamma_{1}=8,\gamma_{2}=\tfrac{91}{2},\gamma_{3}=\tfrac{629}{3},\gamma_{4}=\tfrac{20003}{24},\gamma_{5}=\tfrac{8893}{3},\gamma_{6}=\tfrac{6953959}{720},$
(4.13)

and where

$\displaystyle g_{0}(n)$ $\displaystyle=-(1+n)^{3}(2+n),$ $\displaystyle
g_{1}(n)$ $\displaystyle=(2+n)^{2}(54+40n+7n^{2}),$ $\displaystyle g_{2}(n)$
$\displaystyle=-(2614+3203n+1449n^{2}+287n^{3}+21n^{4}),$ $\displaystyle
g_{3}(n)$ $\displaystyle=10262+10108n+3686n^{2}+590n^{3}+35n^{4},$
$\displaystyle g_{4}(n)$ $\displaystyle=-(4+n)(4430+2688n+535n^{2}+35n^{3}),$
$\displaystyle g_{5}(n)$ $\displaystyle=(4+n)(5+n)(738+251n+21n^{2}),$
$\displaystyle g_{6}(n)$ $\displaystyle=-(4+n)(5+n)(6+n)(48+7n),$
$\displaystyle g_{7}(n)$ $\displaystyle=(n+4)(n+5)(n+6)(n+7).$

Hence, the first $N$ terms in the power series can be enumerated with the help
of a computer. The relative error of the $n^{\text{th}}$ term is shown in
Figure 2 where it can be seen that the convergence of
$\log[E_{n}^{3}]/\sqrt{n}$ to its limiting value is extremely slow.

$200$$300$$400$$500$$600$$700$$800$$900$$1{,}000$$-1.8$$-2$$-2.2$$-2.4$$-2.6$$\sqrt{n}$$\log_{10}\mathopen{}\mathclose{{}\left|\frac{\log
E[X_{n}^{3}]-6.922\sqrt{n}}{6.922\sqrt{n}}}\right|$ Figure 2. The convergence
of $E[X_{n}^{3}]/\sqrt{n}$ to its limiting value from Eq. (4.11).

Performing a linear fit to the terms $(\sqrt{n},\log\gamma_{n})_{900\leq n\leq
1000}$ to compare with [3] returns a slope of roughly $6.83$, which means that
[3] must indeed be an underestimate. Extending our power series to $N=10^{6}$
provides a slope of $6.92$, in agreement with Eq. (4.11).

##  5\. Proofs of the main steps

###  5.1. Proof of Lemma 3.1

We take $n\to n+1$ in Eq. (3.2), then apply Eq. (1.1) raised to the power of
$m$ expanded out with the binomial theorem to obtain

$C_{n+1}^{(m)}(q;\mathfrak{p})=\sum_{\beta=0}^{q}\binom{q}{\beta}\sum_{\alpha_{1}=0}^{n}\sum_{\alpha_{2}=0}^{n}\cdots\sum_{\alpha_{|\mathfrak{p}|}=0}^{n}E\mathopen{}\mathclose{{}\left[X_{n}^{\beta}X_{U(n)}^{q-\beta}\prod_{\ell=1}^{|\mathfrak{p}|}X^{p_{\ell}}_{\alpha_{\ell}}}\right].$
(5.1)

We apply the usual method of conditioning on the value of $U(n)$ inside the
expectation. We expand the result onto three lines which will be analyzed
sequentially;

$\displaystyle C_{n+1}^{(m)}(q;\mathfrak{p})$
$\displaystyle=\sum_{\alpha_{1}=0}^{n}\sum_{\alpha_{2}=0}^{n}\cdots\sum_{\alpha_{|\mathfrak{p}|}=0}^{n}E\mathopen{}\mathclose{{}\left[X_{n}^{q}\prod_{\ell=1}^{|\mathfrak{p}|}X^{p_{\ell}}_{\alpha_{\ell}}}\right]$
(5.2a)
$\displaystyle+\frac{1}{n+1}\sum_{\beta=1}^{q-1}\binom{q}{\beta}\sum_{\alpha_{1}=0}^{n}\sum_{\alpha_{2}=0}^{n}\cdots\sum_{\alpha_{|\mathfrak{p}|}=0}^{n}\sum_{\gamma=0}^{n}E\mathopen{}\mathclose{{}\left[X_{n}^{\beta}X_{\gamma}^{q-\beta}\prod_{\ell=1}^{|\mathfrak{p}|}X^{p_{\ell}}_{\alpha_{\ell}}}\right]$
(5.2b)
$\displaystyle+\frac{1}{n+1}\sum_{\alpha_{1}=0}^{n}\sum_{\alpha_{2}=0}^{n}\cdots\sum_{\alpha_{|\mathfrak{p}|}=0}^{n}\sum_{\gamma=0}^{n}E\mathopen{}\mathclose{{}\left[X_{\gamma}^{q}\prod_{\ell=1}^{|\mathfrak{p}|}X^{p_{\ell}}_{\alpha_{\ell}}}\right].$
(5.2c)

First, we examine the term in Eq. (5.2a). We must write this in terms of sums
which extend to $n-1$ rather than $n$ to compare with Eq. (3.2). To see how
this is done, we note that each sum $\sum_{\alpha_{i}=0}^{n}$ can be expanded
to two kinds of terms: one when $\alpha_{i}=n$ and one for
$\sum_{\alpha_{i}=0}^{n-1}$. Therefore, we can write all the sums as a single
sum over every possible binary vector of length $|\mathfrak{p}|$ such that a
$1$ in the $\ell^{\text{th}}$ entry means that $\alpha_{\ell}=n$ and a $0$
means $\sum_{\alpha_{\ell}=0}^{n-1}$. Consider such a binary vector $v$; the
contribution to the total sum from $v$ is

$E\mathopen{}\mathclose{{}\left[X_{n}^{q}\prod_{p_{\ell}\in\mathfrak{p}}X^{p_{\ell}}_{\alpha_{\ell}}}\right]_{v}=\sum_{\ell:v_{\ell}=0}\sum_{\alpha_{\ell}=0}^{n-1}E\mathopen{}\mathclose{{}\left[X_{n}^{q+\sum_{i:v_{i}=1}p_{i}}\prod_{i:v_{i}=0}X^{p_{\ell}}_{\alpha_{\ell}}}\right]$
(5.3)

The sum over all such $v$ can be equivalently expressed as a sum over every
possible subset of $\mathfrak{p}$ (including the empty set) so that the term
in Eq. (5.2a) can be written as

$\sum_{\mathfrak{b}\subseteq\mathfrak{p}}C_{n}^{(m)}\mathopen{}\mathclose{{}\left(q+\lVert\mathfrak{b}\rVert;\mathfrak{p}-\mathfrak{b}}\right).$
(5.4)

The term in Eq. (5.2b) is similar except for the sum over $\gamma$. By
treating $q-\beta$ as a part, we can write it as

$\sum_{\beta=1}^{q-1}\binom{q}{\beta}\sum_{\mathfrak{b}\subseteq\mathfrak{p}+[q-\beta]}\overline{C}_{n}^{(m)}\mathopen{}\mathclose{{}\left(\beta+\lVert\mathfrak{b}\rVert;\mathfrak{p}+[q-\beta]-\mathfrak{b}}\right),$
(5.5)

with reference to the notation of Eq. (3.3a). The term in Eq. (5.2c) requires
more attention. Temporarily ignoring the factor of $1/(n+1)$, we will write

$\displaystyle D_{n}^{(m)}(q;\mathfrak{p})$
$\displaystyle=\sum_{\alpha_{1}=0}^{n}\sum_{\alpha_{2}=0}^{n}\cdots\sum_{\alpha_{|\mathfrak{p}|}=0}^{n}\sum_{\gamma=0}^{n}E\mathopen{}\mathclose{{}\left[X_{\gamma}^{q}\prod_{\ell=1}^{|\mathfrak{p}|}X^{p_{\ell}}_{\alpha_{\ell}}}\right]$
(5.6a)
$\displaystyle=\sum_{\begin{subarray}{c}\mathfrak{b}\subseteq\mathfrak{p}+[q]\\\
\mathfrak{b}\neq\emptyset\end{subarray}}C_{\ell}^{(m)}\mathopen{}\mathclose{{}\left(\lVert\mathfrak{b}\rVert;\mathfrak{p}+[q]-\mathfrak{b}}\right)+D_{n-1}^{(m)}(q;\mathfrak{p}).$
(5.6b)

In Eq. (5.6b), the $\mathfrak{b}=\emptyset$ term has been separated; note that
$\mathfrak{b}=\emptyset$ corresponds to the case where no $\alpha_{i}$ is
equal to $n$ and hence there is no term of the form $X_{n}^{\ell}$ inside the
expectation so no identification with Eq. (3.2) is possible. But then, Eq.
(5.6b) is a linear recurrence relation for $D_{n}^{(m)}$ which can be directly
solved to obtain

$D_{n}^{(m)}(q;\mathfrak{p})=\sum_{\ell=0}^{n}\sum_{\begin{subarray}{c}\mathfrak{b}\subseteq\mathfrak{p}+[q]\\\
\mathfrak{b}\neq\emptyset\end{subarray}}C_{\ell}^{(m)}\mathopen{}\mathclose{{}\left(\lVert\mathfrak{b}\rVert+[q];\mathfrak{p}-\mathfrak{b}}\right),$
(5.7)

and so, by again referencing the notation of Eq. (3.3a), the term in Eq.
(5.2c) can be written as

$\sum_{\begin{subarray}{c}\mathfrak{b}\subseteq\mathfrak{p}+[q]\\\
\mathfrak{b}\neq\emptyset\end{subarray}}\hat{C}_{\ell}^{(m)}\mathopen{}\mathclose{{}\left(\lVert\mathfrak{b}\rVert;\mathfrak{p}+[q]-\mathfrak{b}}\right).$
(5.8)

Collecting Eq. (5.4), (5.5) and (5.8) gives Eq. (3.4).

###  5.2. Proof of Lemma 3.4

To establish this result, we must investigate the properties of
$\mathbf{C}_{1},\mathbf{C}_{2}$ and $\mathbf{C}_{3}$. The sizes and entries of
these matrices all depend on $m$, and the arguments in this section apply for
any $m\geq 1$. As per Eq. (3.13), the matrix
$(\mathbf{I}-z\mathbf{C}_{1})^{-1}$ will be important, we therefore define

$\mathbf{Q}=\mathbf{I}-z\mathbf{C}_{1},$ (5.9)

note that $\mathbf{Q}$ is lower triangular with diagonal entries $1-z$ by
Proposition. 3.3. To understand $\mathbf{Q}^{-1}$, we will use the following
lemma, which is a standard result from linear algebra for the inverse a
triangular matrix. The proof of Proposition 5.1 is in Appendix A.

######  Proposition 5.1.

Let $X$ be an $n\times n$ lower triangular matrix,

$X=\begin{pmatrix}c_{1,1}&0&0&\dots&0\\\ c_{2,1}&c_{2,2}&0&\dots&0\\\
c_{3,1}&c_{3,2}&c_{3,3}&\dots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
c_{n,1}&c_{n,2}&c_{n,3}&\dots&c_{n,n}\end{pmatrix},$ (5.10)

where $c_{i,i}>0$ for each $1\leq i\leq n$. Then, the entries $a_{i,j}$ of
$X^{-1}$ obey

$a_{i,j}=\begin{cases}\frac{1}{c_{i,i}}\mathopen{}\mathclose{{}\left(\delta_{i,j}-\sum_{k=1}^{i-j}c_{i,i-k}a_{i-k,j}}\right)&1\leq
j\leq i\leq n\\\ 0&\text{otherwise}.\end{cases}$ (5.11)

We will elaborate on the interpretation of Proposition 5.1 in terms of the
indices $(q;\mathfrak{p})$. In particular, consider the sum appearing in Eq.
(5.11),

$\sum_{k=1}^{i-j}c_{i,i-k}a_{i-k,j}.$ (5.12)

The term $c_{i,i-k}$ relates to the dependence of $(q;\mathfrak{p})_{i}$ in
the canonical ordering on $(q;\mathfrak{p})_{i-k}$. The term $a_{i-k,j}$
relates to the dependence of the $(q;\mathfrak{p})_{i-k}$ on
$(q;\mathfrak{p})_{j}$, albeit in in the inverse of the matrix. The sum is
shown schematically for a particular term in the $m=4$ case in Fig. 3.

$(4,\emptyset)$$(3,\\{1\\})$$j=2$$k=5$$(2,\\{2\\})$$k=4$$(1,\\{3\\})$$k=3$$(2,\\{1,1\\})$$k=2$$(1,\\{2,1\\})$$k=1$$(1,\\{1,1,1\\})$$i=7$$c_{i,i-2}$$a_{i-2,j}$
Figure 3. In the $m=4$ case, the terms contributing to Eq. (5.12) for the
calculation of $a_{6,2}$. The sum contains five terms, one for each $1\leq
k\leq 5$. The $k=2$ term is shown in this diagram; the arrows show how the
indices $i,j$ are related to $(q;\mathfrak{p})$ indices in the canonical
ordering.

We are now ready to study the inverse of $\mathbf{Q}$.

######  Proposition 5.2.

For the matrix $\mathbf{Q}$ defined in Eq. (5.9), it holds for all allowed
indices $i,j$ and for $0\leq z<1$ that

  1. (i)

$\mathbf{Q}_{i,j}^{-1}=0$ if and only if $\mathbf{Q}_{i,j}=0$

  2. (ii)

$\mathbf{Q}_{i,j}^{-1}\geq 0$

  3. (iii)

If $\mathbf{Q}_{i,j}^{-1}\neq 0$, then
$\mathbf{Q}_{i,j}^{-1}=(1-z)^{-(|\mathfrak{p}_{i}|-|\mathfrak{p}_{j}|+1)}f_{i,j}(z)$
where $f_{i,j}(z)$ is a polynomial in $z$ such that $f_{i,j}\neq 0$.

###### Proof.

By Proposition 3.3, we can apply the result of Proposition 5.1 to the entries
of $\mathbf{Q}^{-1}$.

We will begin by showing ((i)) by induction on the rows of $\mathbf{Q}^{-1}$.
When there are two rows we have

$\mathbf{Q}^{-1}=\begin{pmatrix}1-z&0\\\
-Q_{1,2}&1-z\end{pmatrix}^{-1}=\begin{pmatrix}(1-z)^{-1}&0\\\
Q_{1,2}(1-z)^{-2}&(1-z)^{-1}\end{pmatrix}.$ (5.13)

which provides the base case for all the assertions of Lemma 3.4 since
$\mathbf{Q}_{1,2}>0$ via the fact that every $\mathfrak{p}$ has the empty set
as a multisubset. Proposition 5.1, shows directly that
$\mathbf{Q}_{i,j}^{-1}=0$ for $j>i$ and $\mathbf{Q}_{i,i}^{-1}=1/(1-z)$, hence
we will restrict our attention to the case $j<i$. On the induction step, we
have

$\mathbf{Q}_{n,j}^{-1}=-\frac{1}{1-z}\sum_{k=1}^{n-j}\mathbf{Q}_{n,n-k}\mathbf{Q}^{-1}_{n-k,j}.$
(5.14)

Suppose first that $\mathbf{Q}_{n,j}=0$; we will examine the $k^{\text{th}}$
term in the sum in Eq. 5.14. If $\mathbf{Q}_{n,n-k}=0$ in Eq. 5.14 then that
term contributes zero to the total sum. If $\mathbf{Q}_{n,n-k}\neq 0$, it
means that $\mathfrak{p}_{n-k}\subset\mathfrak{p}_{n}$. But then, since
$\mathbf{Q}_{n,j}=0$ we have that
$\mathfrak{p}_{j}\not\subset\mathfrak{p}_{n}$ which implies that
$\mathfrak{p}_{j}\not\subset\mathfrak{p}_{n-k}$ and hence that
$\mathbf{Q}_{n-k,j}=0$. By the induction hypothesis, this implies that
$\mathbf{Q}^{-1}_{n-k,j}=0$ and so every term in the sum is zero and we have
$\mathbf{Q}_{n,j}=0\implies\mathbf{Q}^{-1}_{n,j}=0$. Then we can apply Lemma
5.1 with the roles of $\mathbf{Q}$ and $\mathbf{Q}^{-1}$ switched and the same
argument above to obtain $\mathbf{Q}_{n,j}^{-1}=0\implies\mathbf{Q}_{n,j}=0$.
This completes the induction and the proof of ((i)).

To establish ((ii)), we perform a similar induction. Since $\mathbf{C}_{1}\geq
0$, we have $\mathbf{Q}_{i,j}<0$ for $j<i$ and $1\leq z<1$. The induction
hypothesis provides that the sum in Eq. 5.14 contains only non-positive terms
and so $\mathbf{Q}^{-1}_{n,j}$ itself is non-negative.

Finally, we establish ((iii)) by again using similar induction. Since
$\mathbf{Q}_{n,j}^{-1}\neq 0$ by assumption, ((i)) provides that
$\mathfrak{p}_{j}\subset\mathfrak{p}_{n}$. The multiset
$\mathfrak{p}_{n}-\mathfrak{p}_{j}$ contains at least one element, hence by
the requirement of the canonical ordering that generating functions be sorted
by ascending $|\mathfrak{p}|$, there is a $1\leq k\leq n-j$ such that
$\mathfrak{p}_{j}\subseteq\mathfrak{p}_{n-k}\subset\mathfrak{p}_{n}$ and
$|\mathfrak{p}_{n-k}|=|\mathfrak{p}_{n}|-1$. This implies that
$\mathbf{Q}_{n,n-k}>0$ and so by the induction hypothesis, the contribution to
$\mathbf{Q}_{n,j}^{-1}$ in Eq. 5.14 from this term is

$-\frac{1}{1-z}\mathopen{}\mathclose{{}\left[\frac{\mathbf{Q}_{n,n-k}f_{n-k,k}(z)}{(1-z)^{(|\mathfrak{p}_{n}|-1)-|\mathfrak{p}_{j}|+1}}}\right]=\frac{g_{n-k,k}(z)}{(1-z)^{|\mathfrak{p}_{n}|-|\mathfrak{p}_{j}|+1}},$
(5.15)

where $g_{n-k,k}(z)$ is a polynomial in $z$ such that $g_{n-k,k}(1)\neq 0$. By
the same argument, for each multisubset of
$\mathfrak{p}_{n}-\mathfrak{p}_{j}$, there will be a $1\leq k\leq n-j$ where
such that $\mathfrak{p}_{j}\subseteq\mathfrak{p}_{n-k}\subset\mathfrak{p}_{n}$
whose contribution to $\mathbf{Q}_{n,j}^{-1}$ has the same form as Eq. (5.15).
But since $|\mathfrak{p}_{n-k}|\leq|\mathfrak{p}_{n}|-1$, the exponent of
$(1-z)$ is at most $|\mathfrak{p}_{n}|-|\mathfrak{p}_{j}|+1$. Therefore,
$\mathbf{Q}_{n,j}^{-1}$ is of the required form which completes the induction
and the proof. ∎

Referring to Eq. (3.14), Proposition 5.2 implies that we can write

$\mathbf{P}^{-1}(z)\mathbf{Q}^{-1}(z)\mathbf{P}(z)=\frac{1}{1-z}\mathbf{B}_{1}+\mathbf{B}_{2}(z),$
(5.16)

where $\mathbf{B}_{1}$ is a matrix with non-negative entries arranged in a
structure identical to $\mathbf{C}_{1}$ and $\mathbf{B}_{2}(z)$ is a matrix
whose entries are polynomials in $(1-z)$. Define

$\mathbf{D}(z)=\mathbf{P}^{-1}(z)\mathopen{}\mathclose{{}\left[\mathbf{C}_{1}+\mathbf{C}_{2}+\frac{1}{1-z}\mathbf{C}_{3}}\right]\mathbf{P}(z),$
(5.17)

then we can write

$\mathbf{P}^{-1}(z)\mathbf{R}(z)\mathbf{P}(z)=\mathopen{}\mathclose{{}\left(\frac{1}{1-z}\mathbf{B}_{1}+\mathbf{B}_{2}(z)}\right)\mathbf{D}(z)$
(5.18)

We will now study the matrix $\mathbf{D}.$ The effect of a similarity
transformation by $\mathbf{P}(z)$ on a matrix with entries $(a_{i,j})$ is to
transform the entries according to

$a_{i,j}\to a_{i,j}(1-z)^{|\mathfrak{p}_{i}|-|\mathfrak{p}_{j}|}.$ (5.19)

From Eq. (3.4), and by the same arguments as in the proof of Proposition 3.3,
if entry $(i,j)$ in $\mathbf{C}_{2}$ is nonzero, then
$|\mathfrak{p}_{i}|-|\mathfrak{p}_{j}|\geq-1$. Similarly, if entry $(i,j)$ in
$\mathbf{C}_{3}$ is nonzero, then $|\mathfrak{p}_{i}|-|\mathfrak{p}_{j}|\geq
0$. Hence, we can write

$\mathbf{D}(z)=\frac{1}{1-z}\mathbf{E}+\mathbf{F}(z)$ (5.20)

where the entries $\mathbf{F}(z)$ are polynomials in $(1-z)$ and the entries
of $\mathbf{E}$ are obtained by restricting the sum in Eq. (3.4b) to the term
$\mathfrak{b}=\emptyset$ and by restricting the sum in Eq. (3.4c) to terms
where $|\mathfrak{b}|=1$. Explicitly, in a row indexed by $(q;\mathfrak{p})$,
the nonzero entries are either indexed by $(\beta;\mathfrak{p}+[q-\beta])$ for
some $1\leq\beta\leq q-1$ or by $(p;\mathfrak{p}-[p])$ where
$[p]\subseteq\mathfrak{p}+[q]$. Note that

$\mathbf{P}^{-1}(z)\frac{\text{d}}{\text{d}z}\mathbf{P}=\frac{1}{1-z}\text{diag}\\{|\mathfrak{p}_{1}|,|\mathfrak{p}_{2}|,\dots|\mathfrak{p}_{\mathfrak{C}(m)}|\\}\coloneqq\frac{1}{1-z}\mathbf{G},$
(5.21)

hence by referring to Eqs. (3.14), and making the substitution $u=1/(1-z)$ we
can now write

$h^{\prime}(u)=\mathopen{}\mathclose{{}\left[\mathbf{B}_{1}\mathbf{E}+\frac{1}{u}\mathopen{}\mathclose{{}\left(\mathbf{B}_{1}\mathbf{F}(u)+\mathbf{B}_{2}(u)\mathbf{E}-\mathbf{G}}\right)+\frac{1}{u^{2}}\mathbf{B}_{2}(u)\mathbf{F(u)}}\right]h(u),$
(5.22)

where $\mathbf{F}(u)$ and $\mathbf{B}_{2}$ can each be written as a finite sum
of constant matrices multiplied by non-negative powers of $1/u$. Hence, to
complete the proof of Lemma 3.4 we need only establish that
$\mathbf{B}_{1}\mathbf{E}$ is a primitive matrix with the appropriate bounds
on its Perron-Frobenius eigenvalue, which is the subject of the following
proposition. Then, Lemma 3.4 holds with
$\mathbf{M}_{0}=\mathbf{B}_{1}\mathbf{E}$.

######  Proposition 5.3.

The matrix $\mathbf{B}_{1}\mathbf{E}$ is primitive, where $\mathbf{B}_{1}$ is
defined in Eq. (5.16) and $\mathbf{E}$ is defined in Eq. (5.20). The Perron-
Frobenius eigenvalue of $\mathbf{B}_{1}\mathbf{E}$ is bounded above by
$m^{2}2^{m}e^{6\sqrt{m}}$.

###### Proof.

To show that $\mathbf{B}_{1}\mathbf{E}$ is primitive we will show that the
adjacency graph defined by $\mathbf{B}_{1}\mathbf{E}$ with vertices labelled
by the $(q,\mathfrak{p})_{i}$ is strongly connected. Recall that the entries
of $\mathbf{E}$ imply that $(q;\mathfrak{p})$ can reach either
$(\beta;\mathfrak{p}+[q-\beta])$ for some $1\leq\beta\leq q-1$ or by
$(p;\mathfrak{p}-[p])$ where $[p]\subseteq\mathfrak{p}+[q]$ and the entries of
$\mathbf{B}_{1}$ imply that $(q;\mathfrak{p})$ can reach
$(q+\lVert\mathfrak{b}\rVert,\mathfrak{p}-\mathfrak{b})$ for any
$\mathfrak{b}\subseteq\mathfrak{p}$. To demonstrate that the associated
adjacency graph is connected, we will show that there is a bidirectional path
between any $(q;\mathfrak{p})$ and $(m,\emptyset)$. Suppose first that we are
starting at some $(q;\mathfrak{p})$ and want to find a path to
$(m;\emptyset)$. Then, $\mathbf{E}$ provides a path from $(q;\mathfrak{p})$ to
itself and $\mathbf{B}_{1}$ provides a path from $(m;\emptyset)$. Now we would
like to find a path starting from $(m;\emptyset)$ back to $(q;\mathfrak{p})$.
First, $\mathbf{E}$ provides a path from $(m;\emptyset)$ to
$(m-p_{1};[p_{1}])$ then $\mathbf{B}_{1}$ provides a path from
$(m-p_{1};[p_{1}])$ to itself. Iterating this procedure, $\mathbf{E}$ provides
a path from $(m-p_{1},[p_{1}])$ to $(m-p_{1}-p_{2},[p_{1},p_{2}])$ and
$\mathbf{B_{1}}$ provides a path from $(m-p_{1}-p_{2},[p_{1},p_{2}])$ to
itself. It’s clear that that there is a $1\leq k\leq\mathfrak{C}(m)$ such that
$q=m-\sum_{\ell=1}^{k}p_{\ell}$ and hence there will be a path from
$(m;\emptyset)$ to $(q;\mathfrak{p})$ as required. This completes the proof
that the matrix $\mathbf{B}_{1}\mathbf{E}$ is primitive.

Now we will establish the appropriate bounds on the Perron-Frobenius
eigenvalue $\lambda$ of $\mathbf{B}_{1}\mathbf{E}$. To do this, we will use
the fact that $\lambda$ is bounded above by the maximum sum of the columns of
$\mathbf{B}_{1}\mathbf{E}$. By Lemma 3.1, the largest entry that can appear in
$\mathbf{C}_{1}$ is $\binom{m}{\lfloor m/2\rfloor}\leq 2^{m}$. Since
$\mathfrak{C}(m)\leq e^{3\sqrt{m}}$ [12], Proposition 5.1, implies via
induction that the largest entry of $\mathbf{B}_{1}$ is bounded above by
$2^{me^{3\sqrt{m}}}$. Similarly, the largest entry that can appear in
$\mathbf{E}$ is $\binom{m}{\lfloor m/2\rfloor}\leq 2^{m}$. It follows that
each entry of $\mathbf{B}_{1}\mathbf{E}$ is bounded above by
$2^{(m+1)e^{3\sqrt{m}}}e^{3\sqrt{m}}$ and hence that $\lambda$ is bounded
above by
$2^{\mathopen{}\mathclose{{}\left(\tfrac{m}{2}+1}\right)e^{3\sqrt{m}}}$. ∎

###  5.3. Proof of Lemma 3.5

To prove Lemma 3.5, we will apply two theorems from the theory of
perturbations of linear systems. The following theorem is originally due to
Hartman and Wintner [13] and provides to us the required asymptotic form of
the solutions to our differential equation system. We quote the variation
appearing in [19].

######  Theorem 5.4 ([19], Theorem 1).

Let $\mathbf{A}$ be a constant $d\times d$ matrix and let $y(t)=(y_{i})_{1\leq
i\leq d}$ be a column vector. Consider the system of differential equations

$y^{\prime}(t)=\mathopen{}\mathclose{{}\left(\mathbf{A}+\mathbf{B}(t)}\right)y(t),$
(5.23)

where $\mathbf{B}(t)$ is a continuous $d\times d$ matrix for $t\geq 0$. Let
$\lambda$ be a simple eigenvalue of $\mathbf{A}$ such that no other eigenvalue
of $\mathbf{A}$ has the same real part as $\lambda$ and let $v$ be the
eigenvector associated to this eigenvalue. If $B(t)\in L^{2}(t_{0},\infty)$
then, for $t_{0}$ large enough, Eq. (5.23) has a solution on $[t_{0},\infty)$
such that

$y(t)=\exp\mathopen{}\mathclose{{}\left[\mu(t-t_{0})+\int_{t_{0}}^{t}\delta(\tau)\,\text{d}\tau}\right]\mathopen{}\mathclose{{}\left[v+o(1)}\right]\quad\text{as
}t\to\infty,$ (5.24)

where

$\delta(t)=\frac{(\mathbf{B}(t)v)\cdot v}{v\cdot v}.$ (5.25)

The following is Theorem 1.10.1 of [9] modified to our simpler case.

######  Theorem 5.5 ([9], Theorem 1.10.1).

Let $\mathbf{J}=\kappa\mathbf{I}+\mu\mathbf{E}$ be an $n\times n$ matrix of
Jordan type and let $\mathbf{R}(t)$ be an $n\times n$ matrix continuous on
$[t_{0},\infty)$. Consider the system

$y^{\prime}(t)=(J+R(t))y(t).$ (5.26)

Then, if $\mathbf{R}(t)\in L^{2}(t_{0},\infty)$, Eq. (5.26) has solutions
$w_{1}(t),w_{2}(t)\dots w_{n}(t)$ such that for some constants $\alpha_{i}$,

$w_{i}(t)\sim t^{\alpha_{i}}e^{\kappa t}.$ (5.27)

We can now prove Lemma 3.5.

###### Proof of Lemma 3.5.

We will use the fact that the Perron-Frobenius eigenvalue of a matrix is real,
simple and larger than any other eigenvalue in absolute value. We apply
Theorem 5.4 with $\mathbf{A}=\mathbf{M}_{0}$ and
$\mathbf{B}(t)=\sum_{\ell=1}^{n}M_{\ell}t^{-\ell}$. Since $\mathbf{B}(t)$ is
at most of order $1/t$ for $t\geq 1$, all the conditions of Theorem 5.4 are
satisfied, and we obtain that Eq. (3.17) has a solution of the form in Eq.
(3.18).

Next we show that there are no solutions which grow faster than this solution
in the log-asymptotic sense. Transform the system Eq. (3.17) such that
$\mathbf{M}_{0}$ is in Jordan normal form. Then, each Jordan block defines a
subsystem of equations such that the result of Theorem 5.5 applies. Solutions
to these subsystems grow log-asymptotically as fast as the eigenvalue
associated to that Jordan block, the largest of which is the Perron-Frobenius
eigenvalue of $\mathbf{M}_{0}$. ∎

##  6\. Open Problems

Following the computation of the matrix $\mathbf{B}_{1}\mathbf{E}$ in Section
5.2, the matrix $\mathbf{M}_{0}$ of Lemma 3.5 can be computed exactly for
small values of $m$ with the help of a computer. For example, with $m=4$ we
have

$\mathbf{M}_{0}(m=4)=\mathopen{}\mathclose{{}\left(\begin{array}[]{ccccccc}1&4&6&4&0&0&0\\\
1&5&6&5&3&3&0\\\ 1&4&8&4&0&2&0\\\ 1&5&6&5&0&0&0\\\ 2&10&12&10&7&8&2\\\
2&9&14&9&4&7&0\\\ 6&30&36&30&21&24&10\\\ \end{array}}\right).$ (6.1)

By computing the Perron-Frobenius eigenvalues of the $\mathbf{M}_{0}$ for each
$m$, we can calculate the proportionality constants $c_{m}$ in $\log
E[X_{n}^{m}]\sim c_{m}\sqrt{n}$ exactly, the results are shown in Figure 4.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$0$$10$$20$$30$$40$$50$$60$$70$$m$$c_{m}$
Figure 4. The exact values of $c_{m}$ for which $\log E[X_{n}^{m}]\sim
c_{m}\sqrt{n}$ for small $m$.

Comparing the values of $c_{m}$ in Figure 4 to the bounds provided in Theorem
1.3, it’s clear that the bounds in Theorem 1.3 can be tightened significantly.
Based on numerical evidence, we propose the following conjecture.

######  Conjecture 6.1.

Theorem 1.3 holds with the bounds on $c_{m}$ replaced by the stronger bounds
$\log c_{m}\in[0.3m,0.4m+1]$.

Estimating the size of the Perron-Frobenius eigenvalue based on the bounds
provided by row and column sums does not appear to be strong enough to obtain
this. As observed in [16], the generalization of the Ulam-Kac adder to Eq.
(1.4) still admits first and second moments which grow log-asymptotically as
$\sqrt{n}$. We therefore propose the following additional conjecture, which
can likely be decided with only minor modifications to the arguments in this
paper,

######  Conjecture 6.2.

Theorem 1.3 holds for Eq. (1.4) with the bounds on $c_{m}$ replaced by
$[f_{1}(m,\gamma),f_{2}(m,\gamma)]$ for suitable functions $f_{1},f_{2}$.

Conjecture 1.2 as well as questions related to the first passage behavior of
the Ulam-Kac adder remain completely open. Many history-dependent random
sequences which have been studied recently share the general feature of
defining recurrence relations which contain terms of the form $X_{U(n)}$. We
considered the case where $U(n)$ is a discrete uniform distribution, but
generalizations to non-uniform distributions such as those mentioned in [16]
require further investigation.

## Acknowledgements

The author would like to thank Benedek Valkó for providing several very
helpful comments on an earlier draft of this paper as well as Yu-Qiu Zhao for
providing a number of references on the asymptotic theory of difference
equations.

##  Appendix A Supplementary Material

###### Proof of Lemma 5.1.

The proof is by induction on $n$. For the base case $n=1$, we have
$X^{-1}=(1/c_{1,1})$ and Eq. (5.11) reads $a_{1,1}=1/c_{1,1}$. The induction
step proceeds by Gaussian elimination of the $n^{\text{th}}$ row with the
first $n-1$ rows already reduced. We then have

$\mathopen{}\mathclose{{}\left[\begin{array}[]{ccccc|ccccc}1&0&0&\dots&0&a_{1,1}&0&0&\dots&0\\\
0&1&0&\dots&0&a_{2,1}&a_{2,2}&0&\dots&0\\\
0&0&1&\dots&0&a_{3,1}&a_{3,2}&a_{3,3}&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
c_{n,1}&c_{n,2}&c_{n,3}&\dots&c_{n,n}&0&0&0&\dots&1\end{array}}\right]$ (A.1)

We see that to reduce $c_{n,j}$ for $1\leq j\leq n-1$, we multiply the
$j^{\text{th}}$ row by $-c_{n,j}$ and add it to the $n^{\text{th}}$ row. This
contributes a term $-c_{n,j}a_{j,k}$ to $a_{n,k}$ for each $1\leq k\leq j$.
Then, to reduce $c_{n,n}$, we divide the $n^{\text{th}}$ row by $c_{n,n}$.
Therefore,

$a_{n,j}=\frac{1}{c_{n,n}}\mathopen{}\mathclose{{}\left(\delta_{n,j}-\sum_{k=1}^{n-j}c_{n,n-k}a_{n-k,j}}\right),\quad
1\leq j\leq n,$ (A.2)

which completes the induction.

∎

## References

  * [1] H. M. Bahig, Star reduction among minimal length addition chains, Computing, 91 (2011), pp. 335–352.
  * [2] H. M. Bahig and Y. Kotb, An efficient multicore algorithm for minimal length addition chains, Computers, 8 (2019), p. 23.
  * [3] E. Ben-Naim and P. L. Krapivsky, Growth and structure of stochastic sequences, Journal of Physics A: Mathematical and General, 35 (2002), p. L557–L563.
  * [4] C. M. Bender, S. Orszag, and S. A. Orszag, Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory, vol. 1, Springer Science & Business Media, 1999.
  * [5] J. Benito, N. C. Karayiannis, and M. Laso, Confined polymers as self-avoiding random walks on restricted lattices, Polymers, 10 (2018), p. 1394.
  * [6] A. Brauer, On addition chains, Bulletin of the American mathematical Society, 45 (1939), pp. 736–739.
  * [7] P. Clifford and D. Stirzaker, History-dependent random processes, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 464 (2008), pp. 1105–1124.
  * [8] C. Domb, Self-Avoiding Walks on Lattices, John Wiley & Sons, Ltd, 1969, pp. 229–259.
  * [9] M. S. P. Eastham et al., The asymptotic solution of linear differential systems: applications of the Levinson theorem, no. 4, Oxford University Press, 1989.
  * [10] P. Flajolet and R. Sedgewick, Analytic combinatorics, cambridge University press, 2009.
  * [11] K. F. Freed, Polymers as self-avoiding walks, The Annals of Probability, (1981), pp. 537–554.
  * [12] H. Gupta, On an asymptotic formula in partitions, in Proceedings of the Indian Academy of Sciences-Section A, vol. 16, Springer India, 1942, pp. 101–102.
  * [13] P. Hartman and A. Wintner, Asymptotic integrations of linear differential equations, American Journal of Mathematics, 77 (1955), pp. 45–86.
  * [14] M. Kac, A history-dependent random sequence defined by ulam, Advances in Applied Mathematics, 10 (1989), pp. 270 – 277.
  * [15] M. Kadir, M. Mohamed, R. Mohamad, M. Mamat, and A. Muhammed, Performance comparison of some addition chain methods based on integer family, in International Conference on Information Science and Applications, Springer, 2018, pp. 211–218.
  * [16] I. Krasikov, G. Rodgers, and C. Tripp, Growing random sequences, Journal of Physics A General Physics, 37 (2004).
  * [17] D. MacDonald, S. Joseph, D. Hunter, L. Moseley, N. Jan, and A. Guttmann, Self-avoiding walks on the simple cubic lattice, Journal of Physics A: Mathematical and General, 33 (2000), p. 5973.
  * [18] O. Parreño, P. M. Ramos, N. C. Karayiannis, and M. Laso, Self-avoiding random walks as a model to study athermal linear polymers under extreme plate confinement, Polymers, 12 (2020), p. 799.
  * [19] M. Pituk, The hartman–wintner theorem for functional differential equations, journal of differential equations, 155 (1999), pp. 1–16.
  * [20] , More on poincaré’s and perron’s theorems for difference equations, The Journal of Difference Equations and Applications, 8 (2002), pp. 201–216.
  * [21] A. Roitershtein and Z. Zhou, Distribution tails of a history-dependent random linear recursion, Stochastic Models, 38 (2022), pp. 250–267.
  * [22] A. Scholz, Aufgabe 253, Jahresbericht der deutchen Mathematiker-Vereinigung, 47 (1937), pp. 41–42.
  * [23] A. Schönhage, A lower bound for the length of addition chains, Theoretical Computer Science, 1 (1975), pp. 1–12.
  * [24] K. B. Stolarsky, A lower bound for the scholz-brauer problem, Canadian Journal of Mathematics, 21 (1969), pp. 675–683.
  * [25] G. Szeg, Orthogonal polynomials, vol. 23, American Mathematical Soc., 1939.
  * [26] S. M. Ulam et al., Analogies between analogies: the mathematical reports of SM Ulam and his Los Alamos collaborators, vol. 10, Univ of California Press, 1990.
  * [27] J. Wimp and D. Zeilberger, Resurrecting the asymptotics of linear recurrences, Journal of mathematical analysis and applications, 111 (1985), pp. 162–176.