text:
\begin{document}

\maketitle

\begin{abstract}

A recent Graph Neural Network (GNN) approach for learning to branch has been shown to successfully reduce the running time of branch-and-bound (\bnb{}) algorithms for Mixed Integer Linear Programming (MILP).
While the GNN relies on a GPU for inference, MILP solvers are purely CPU-based.
This severely limits its application as many practitioners may not have access to high-end GPUs.
In this work, we ask two key questions.
First, in a more realistic setting where only a CPU is available, is the GNN model still competitive?
Second, can we devise an alternate computationally inexpensive model that retains the predictive power of the GNN architecture?
We answer the first question in the negative, and address the second question by proposing a new hybrid architecture for efficient branching on CPU machines.
The proposed architecture combines the expressive power of GNNs with computationally inexpensive multi-layer perceptrons (MLP) for branching.
We evaluate our methods on four classes of MILP problems, and show that they lead to up to 26\% reduction in solver running time compared to state-of-the-art methods without a GPU, while extrapolating to harder problems than it was trained on. The code for this project is publicly available at \url{https://github.com/pg2455/Hybrid-learn2branch}.
\end{abstract}

\section{Introduction}

Mixed-Integer Linear Programs (MILPs) arise naturally in many decision-making problems such as auction design~\citep{abrache2007combinatorial}, warehouse planning~\citep{elson1972site}, capital budgeting~\citep{capbudgeting} or scheduling~\citep{floudas2005mixed}.
Apart from a linear objective function and linear constraints, some decision variables of a MILP are required to take integral values, which makes the problem NP-hard~\citep{book/PapadimitriouS82}.

Modern mathematical solvers typically employ the \bnb{} algorithm~\citep{journal/econometrica/LandD60} to solve general MILPs to global optimality.
While the worst-case time complexity of \bnb{} is exponential in the size of the problem ~\citep{books/Wolsey98}, it has proven efficient in practice, leading to wide adoption in various industries.
At a high level, \bnb{} adopts a divide-and-conquer approach that consists in recursively partitioning the original problem into a tree of smaller sub-problems, and solving linear relaxations of the sub-problems until an integral solution is found and proven optimal.

Despite its apparent simplicity, there are many practical aspects that must be considered for \bnb{} to perform well~\citep{achterberg_thesis}; such decisions will affect the search tree, and ultimately the overall running time.
These include several decision problems~\citep{journals/top/LodiZ2017} that arise during the execution of the algorithm, such as \textit{node selection}: which sub-problem do we analyze next?; and \textit{variable selection} (a.k.a. branching): which decision variable must be used (branched on) to partition the current sub-problem?
While such decisions are typically made using hard-coded expert heuristics which are implemented in modern solvers, more and more attention is given to statistical learning approaches for replacing and improving upon those heuristics \citep{conf/nips/HeDE14,Alvarez2014ASM, confs/aaai/KhalilBND16,gasse2019exact,zarpellon2020parameterizing}.
An extensive review of different approaches at the intersection of statistical learning and combinatorial optimization is given in \citet{arxiv/BengioLP18}.

\input{inputs/times_arch}
Recently, \citet{gasse2019exact} proposed to tackle the variable selection problem in \bnb{} using a Graph Neural Network (GNN) model.
The GNN exploits the bipartite graph formulation of MILPs together with a shared parametric representation, thus allowing it to model problems of arbitrary size.
Using imitation learning, the model is trained to approximate a very good but computationally expensive ``expert" heuristic named \textit{strong branching} \citep{techreport/ApplegateBCC95}.
The resulting branching strategy is shown to improve upon previously proposed approaches for branching on several MILP problem benchmarks, and is competitive with state-of-the-art \bnb{} solvers.

We note that one limitation of this approach, with respect to general \bnb{} heuristics, is that the resulting strategy is only tailored to the class of MILP problems it is trained on.
This is very reasonable in our view, as practitioners usually only care about solving very specific problem types at any time.

While the GNN model seems particularly suited for learning branching strategies, one drawback is a high computational cost for inference, i.e., choosing the branching variable at each node of the \bnb{} tree.
In~\citet{gasse2019exact}, the authors use a high-end GPU card to speed-up the GNN inference time, which is a common practice in deep learning but is somewhat unrealistic for MILP practitioners.
Indeed, commercial MILP solvers rely solely on CPUs for computation, and the GNN model from \citet{gasse2019exact} is not competitive on CPU-only machines, as illustrated in Figure~\ref{fig:times_arch}.
There is indeed a trade-off between the quality of the branching decisions made and the time spent obtaining those decisions.
This trade-off is well-known in MILP community~\citep{achterberg_thesis}, and has given rise to carefully balanced strategies designed by MILP experts, such as \textit{hybrid branching}~\citep{confs/cpaior/AchterbergB09}, which derive from the computationally expensive \textit{strong branching} heuristic~\citep{techreport/ApplegateBCC95}.

In this paper, we study the time-accuracy trade-off in learning to branch with the aim of devising a model that is both computationally inexpensive and accurate for branching.
To this end, we propose a hybrid architecture that uses a GNN model only at the root node of the \bnb{} tree and a weak but fast predictor, such as a simple Multi-Layer Perceptron (MLP), at the remaining nodes.
In doing so, the weak model is enhanced by high-level structural information extracted at the root node by the GNN model.
In addition to this new hybrid architecture, we experiment and evaluate the impact of several variants to our training protocol for learning to branch, including: (i) end-to-end training~\citep{glasmachers2017limits, bojarski2016end}, (ii) knowledge distillation~\citep{hinton2015distilling}, (iii) auxiliary tasks~\citep{liebel2018auxiliary}, and (iv) depth-dependent weighting of the training loss for learning to branch, an idea originally proposed by  \citet{conf/nips/HeDE14} in the context of node selection.

We evaluate our approach on large-scale MILP instances from four problem families: Capacitated Facility Location, Combinatorial Auctions, Set Covering, and Independent Set.
We demonstrate empirically that our combination of hybrid architecture and training protocol results in state-of-the-art performance in the realistic setting of a CPU-restricted machine.
While we observe a slight decrease in the predictive performance of our model with respect to the original GNN from~\citep{gasse2019exact}, its reduced computational cost still allows for a reduction of up to 26\% in overall solving time on all the evaluated benchmarks compared to the default branching strategy of the modern open-source solver SCIP~\citep{GleixnerEtal2018OO}.
Also, our hybrid model preserves the ability to extrapolate to harder problems than trained on, as the original GNN model.

\section{Related Work}

Finding a branching strategy that results in the smallest \bnb{} tree for a MILP is at least as hard--and possibly much harder--than solving the MILP with any strategy.
Still, small trees can be obtained by using a computationally expensive heuristic named \textit{strong branching} (SB)~\citep{techreport/ApplegateBCC95, lindSB}.
The majority of the research efforts in variable selection are thus aimed at matching the performance of SB through faster approximations, via cleverly handcrafted heuristics such as \textit{reliability pseudocost branching}~\citep{ACHTERBERG200542}, and recently via machine learning\footnote{Interestingly, early works on ML methods for branching can be traced back to 2000 \citep[Acknowledgements]{achterberg_thesis}.}~\citep{Alvarez2014ASM,confs/aaai/KhalilBND16, gasse2019exact}. We refer to~\cite{journals/top/LodiZ2017} for an extensive survey of the topic.

\citet{Alvarez2014ASM} and \citet{confs/aaai/KhalilBND16} showed that a fast discriminative classifier such as extremely randomized trees~\citep{geurts2006extremely} or support vector machines~\citep{svmHearst} on hand-designed features can be used to mimic SB decisions.
Subsequently, \citet{gasse2019exact} and \citet{zarpellon2020parameterizing} showed the importance of representation learning for branching.
Our approach, in some sense, combines the superior representation framework of \citet{gasse2019exact} with the computationally cheaper framework of \citet{confs/aaai/KhalilBND16}.
Such hybrid architectures have been successfully used in ML problems such as visual reasoning~\citep{perez2018film}, style-transfer~\citep{dumoulin2016learned}, natural language processing~\citep{srivastava2015highway, dhingra2016gated}, and speech recognition~\citep{kim2017dynamic}.

\section{Preliminaries}

Throughout this paper, we use boldface for vectors and matrices. A MILP is a mathematical optimization problem that combines a linear objective function, a set of linear constraints, and a mix of continuous and integral decision variables. It can be written as:
\begin{equation*}
    \argmin_{\mathbf{x}} \mathbf{c}^{\intercal}\mathbf{x}, \quad \text{s.t.} \quad  \mathbf{A}\mathbf{x} \leq \mathbf{b}, \quad  \mathbf{x} \in \mathbb{Z}^{p} \times \mathbb{R}^{n-p},

\end{equation*}
where $\mathbf{c} \in \mathbb{R}^n$ denotes the cost vector, $\mathbf{A} \in \mathbb{R}^{m\times n}$ the matrix of constraint coefficients, $\mathbf{b} \in \mathbb{R}^m$ is the vector of constant terms of the constraints, and there are $p$ integer variables, $1\leq p \leq n$ .

The \bnb{} algorithm can be described as follows.
One first solves the linear program (LP) relaxation of the MILP, obtained by disregarding the integrality constraints on the decision variables.
If the LP solution $\mathbf{x}^\star$ satisfies the MILP integrality constraints, or is worse than a known integral solution, then there is no need to proceed further.
If not, then one divides the MILP into two sub-MILPs. This is typically done by picking an integral decision variable that has a fractional value, $i \in \mathcal{C}=\{i \mid x^\star_i \not\in \mathbb{Z}, i \leq p\}$, and create two sub-MILPs with additional constraints $x_i \leq \lfloor x^\star_i \rfloor$ and $x_i \geq \lceil x^\star_i \rceil$, respectively.
The decision variable $i$ that is used to partition the feasible region is called the \textit{branching variable}, while $\mathcal{C}$ denotes the \textit{branching candidates}.
The second step is to select one of the leaves of the tree, and repeat the above steps until all leaves have been processed\footnote{For a more involved description of \bnb{}, the reader is referred to \citet{achterberg_thesis}.}.

In this work, we refer to the first node processed by \bnb{} as the \textit{root node}, which contains the original MILP, and all subsequent nodes containing a local MILP as \textit{tree nodes}, whenever the distinction is required. Otherwise we refer to them simply as nodes.

\section{Methodology}
\label{sec:arch&train}

As mentioned earlier, computationally heavy GNNs can be prohibitively slow when used for branching on CPU-only machines.
In this section we describe our hybrid alternative, which combines the superior inductive bias of a GNN at the root node with a computationally inexpensive model at the tree nodes. We also discuss various enhancements to the training protocol, in order to enhance the performance of the learned models.

\subsection{Hybrid architecture}

A variable selection strategy in \bnb{} can be seen as a scoring function $f$ that outputs a score $s_i \in \mathbb{R}$ for every branching candidate. As such, $f$ can be modeled as a parametric function, learned by ML.
Branching then simply involves selecting the highest-scoring candidate according to $f$:
\begin{equation*}
    \qquad i^\star_f = \argmax_{i \in \mathcal{C}} \mathbf{s}_i
\end{equation*}

We consider two forms of node representations for machine learning models: (i) a graph representation $\mathbf{G} \in \mathcal{G}$, such as the variable-constraint bipartite graph of \citet{gasse2019exact}, where $\mathbf{G} = (\mathbf{V}, \mathbf{E}, \mathbf{C})$, with $\mathbf{V} \in \mathbb{R}^{n \times d_1}$ variable features, $\mathbf{E} \in \mathbb{R}^{n \times m \times d_2}$ edge features, and $\mathbf{C} \in \mathbb{R}^{m \times d_3}$ constraint features; and (ii) branching candidate features
$\mathbf{X} \in \mathbb{R}^{|\mathcal{C}| \times d_4}$, such as those from \citet{confs/aaai/KhalilBND16}, which is cheaper to extract than the first representation.
For convenience, we denote by $\mathcal{X}$ the generic space of the branching candidate features, and by $\mathbf{G}^0$ the graph representation of the root node.
The various features $d_i$ used in this work are detailed in the supplementary materials.

In \bnb{}, structural information in the tree nodes, $\mathbf{G}$, shares a lot of similarity with that of the root node, $\mathbf{G}^0$. Extracting, but also processing that information at every node is an expensive task, which we will try to circumvent. The main idea of our hybrid approach is then to succinctly extract the relevant structural information only once, at the root node, with a parametric model $\text{GNN}(\mathbf{G}^0; \pmb{\theta})$. We then combine in the tree nodes this preprocessed structural information with the cheap candidate features $\mathbf{X}$, using a hybrid model $f := \text{MLP}_\text{HYBRID}(\text{GNN}(\mathbf{G}^0; \pmb{\theta}), \mathbf{X}; \pmb{\phi})$. By doing so, we hope that the resulting model will approach the performance of an expensive but powerful $f := \text{GNN}(\mathbf{G})$, at almost the same cost as an inexpensive but less powerful $f := \text{MLP}(\mathbf{X})$. Figure~\ref{fig:arch-data-type} illustrates the differences between those approaches, in terms of data extraction.

\input{inputs/arch-data-types}

For an exhaustive coverage of the computational spectrum of hybrid models, we consider four ways to enrich the feature space of an MLP via a GNN's output, sumarrized in Table~\ref{tabl:arch-types}.
In CONCAT, we concatenate the candidate's \textit{root representations} $\pmb{\Psi}$ with the features $\mathbf{X}$ at a node.
In FiLM~\citep{perez2018film}, we generate film parameters $\pmb{\gamma}$, $\pmb{\beta}$ from the GNN, for each candidate, which are further used to modulate the hidden layers of the MLP. In details, if $\pmb{h}$ is the intermediate representation of the MLP, it gets linearly modulated as $\pmb{h} \leftarrow \pmb{\beta} \cdot  \pmb{h} + \pmb{\gamma}$.
While both the above architectures have similar computational complexity, it has been shown that FiLM subsumes the CONCAT architecture~\cite{dumoulin2018feature}.
On the other the end of the spectrum lie the most inexpensive hybrid architectures, HyperSVM and HyperSVM-FiLM.
HyperSVM is inspired by HyperNetworks~\citep{ha2016hypernetworks}, and simply consists in a multi-class Support Vector Machine (SVM), whose parameters are predicted by the root GNN.
We chose a simple linear disciminator, for a minimal computational cost.
Finally, in HyperSVM-FiLM we increase the expressivity of HyperSVM with the help of modulations, similar to that of FiLM.
\input{inputs/arch-types}

\subsection{Training Protocol}
We use \textit{strong branching} decisions as ground-truth labels for imitation learning, and collect observations of the form $(\mathbf{G}^0, \mathbf{G}, \mathbf{X}, i^\star_{SB})$.
Thus, the data used for training the model is $\mathcal{D} = \{(\mathbf{G}^0_k, \mathbf{G}_k, \mathbf{X}_k, i^\star_{SB,k}), k=1,2,...,N\}$.
We treat the problem of identifying $i^\star_{SB,k}$ as a classification problem, such that $i^* = \argmax_{i \in \mathcal{C}} f(\mathbf{G}^0, \mathbf{X})$ is the target outcome.
Considering $\mathcal{D}$ as our ground-truth, our objective (\ref{eq:loss-fn}) is to minimize the cross-entropy loss $l$ between $f(\mathbf{G}^0, \mathbf{X}) \in \mathcal{R}^{|\mathcal{C}|}$ and
a one-hot vector with one at the target:
\begin{equation}
   \mathcal{L}(\mathcal{D}; \pmb{\theta}, \pmb{\phi}) = \frac{1}{N}\sum_{k=1}^N l(f(\mathbf{G}^0_k, \mathbf{X}_k ; \pmb{\theta}, \pmb{\phi}), i^\star_{SB,k}).
    \label{eq:loss-fn}
\end{equation}

Performance on unseen instances, or generalization, is of utmost importance when the trained models are used on bigger instances.
The choices in training the aforementioned architectures influence this ability.
In this section, we discuss four such important choices that lead to better generalization.

\subsubsection{End-to-end Training (e2e)}
A $\text{GNN}_{\pmb{\hat{\theta}}}$, with pre-trained parameters $\pmb{\hat{\theta}}$, is obtained using the procedure described in \citet{gasse2019exact}.
We use this pre-trained GNN to extract variable representations at the root node and use it as an input to the MLPs at a tree node (pre).
This results in a considerable performance boost over plain "salt-of-the-earth" MLP used at all tree nodes.
However, going a step further, an end-to-end (e2e) approach involves training the GNN and the MLP together by backpropagating the gradients from a tree node to the root node.
In doing so, e2e training aligns the variable representations at the root node with the prediction task at a tree node.
At the same time, it is not obvious that it should result in a stable learning behavior because the parameters for GNN need to adapt to various tree nodes.
Our experiments explore both pre-training (pre) and end-to-end training (e2e), namely:
\begin{equation}\label{eq:e2e}
    \text{(pre)} \qquad \pmb{\phi}^* = \argmin_{\pmb{\phi}}\mathcal{L}(\mathcal{D}; \pmb{\hat{\theta}}, \pmb{\phi} )
    \text{,} \qquad
    \text{(e2e)} \qquad \pmb{\phi}^*, \pmb{\theta}^* = \argmin_{\pmb{\theta}, \pmb{\phi}}\mathcal{L(\mathcal{D}; \pmb{\theta}, \pmb{\phi})}.
\end{equation}

\subsubsection{Knowledge Distillation (KD)}
Using the outputs of a pre-trained expert model as a soft-target for training a smaller model has been successfully used in model compression~\citep{hinton2015distilling}.
In this way, one aims to learn an inexpensive model that has the same generalization power as an expert.
Thus, for better generalization, instead of training our hybrid architectures with cross-entropy~\citep{Goodfellow-et-al-2016} on ground-truth hard-labels, we study the effect of training with KL Divergence~\citep{kullback1951information} between the outputs of a pre-trained GNN and a hybrid model, namely:
\begin{equation}\label{eq:kd}
   \text{(KD)} \qquad \mathcal{L}_{KD}(\mathcal{D}; \pmb{\theta}, \pmb{\phi}) = \frac{1}{N}\sum_{k=1}^N \text{KL}(f(\mathbf{G}^0_k, \mathbf{X}_k ; \pmb{\theta}, \pmb{\phi}), \text{GNN}_{\pmb{\hat{\theta}}}(\mathbf{G_k})).
\end{equation}

\subsubsection{Auxiliary Tasks (AT)}
An inductive bias, such as GNN, encodes a prior on the way to process raw input.
Auxiliary tasks, on the other hand, inject priors in the model through additional learning objectives, which are not directly linked to the main task.
These tasks are neither related to the final output nor do they require additional training data.
One such auxiliary task is to maximize the diversity in variable representations.
The intuition is that very similar representations lead to very close MLP score predictions, which is not useful for branching.

We minimize a pairwise loss function that ensures maximum separation between the variable representations projected on a unit hypersphere.
We consider two types of objectives for this: (i) Euclidean Distance (ED), and (ii) Minimum Hyperspherical Energy (MHE)~\citep{liu2018learning}, inspired from the well-known Thomson problem~\cite{thomson} in Physics.
While ED separates the representations in the Euclidean space on the hypersphere, MHE ensures uniform distribution over the hypersphere.
Denoting $\hat{\pmb{\psi}_i}$ as the variable representation for the variable $i$ projected on a unit hypersphere and $e_{ij} = ||{\hat{\pmb{\psi}_i} - \hat{\pmb{\psi}_j}}||_2$ as the Euclidean distance between the representations for the variables $i$ and $j$, our new objective function is given as $\mathcal{L}_{AT}(\mathcal{D};\pmb{\theta}, \pmb{\phi}) = \mathcal{L}(\mathcal{D};\pmb{\theta}, \pmb{\phi}) + g(\mathbf{\Psi};\pmb{\theta})$, where
\begin{equation}\label{eq:at}
    \text{(ED)} \qquad  g(\mathbf{\Psi};\pmb{\theta}) = \frac{1}{N^2}\sum_{i,j=1}^{N} e_{ij}^2
    \text{,} \qquad
    \text{(MHE)} \qquad g(\mathbf{\Psi};\pmb{\theta}) = \frac{1}{N^2}\sum_{i,j=1}^{N}\frac{1}{e_{ij}}.
\end{equation}

\subsubsection{Loss Weighting Scheme}
\label{sec:weighting}
The problem of distribution shift is unavoidable in a sequential process like B\&B.
A suboptimal branching decision at a node closer to the root node can have worse impact on the size of the \bnb{} tree as compared to when such a decision is made farther from it.
In such situations, one can use depth (possibly normalized) as a feature, but the generalization on bigger instances is a bit unpredictable as the distribution of this feature might be very different from that observed in the training set.
Thus, we experimented with different depth-dependent formulations for weighting the loss at any node.
Denoting $z_i$ as the depth of a tree node $i$ relative to the depth of the tree, we weight the loss at different tree nodes by $w(z_i)$, making our objective function as
\begin{equation}\label{eq:weights}
   \mathcal{L}(\mathcal{D}; \pmb{\theta}, \pmb{\phi}) = \frac{1}{N}\sum_{k=1}^N w(z_k) \cdot  l(f(\mathbf{G}^0_k, \mathbf{X}_k ; \pmb{\theta}, \pmb{\phi}), i^\star_{SB,k}).
\end{equation}

Specifically, we considered 5 different weighting functions such that all of them have the same end points, i.e., $w(0) = 1.0$ at the root node and $w(1) = e^{-0.5}$ at the deepest node.
Different functions were chosen depending on their intermediate behaviour in between these two points.
We experimented with exponential, linear, quadratic and sigmoidal decay behavior of these functions.
Table~\ref{tab:weights} lists various functions and their mathematical forms considered in our experiments.

\section{Experiments}
\label{sec:experiments}
We follow the experimental setup of \citet{gasse2019exact}, and evaluate each branching strategy across four different problem classes, namely Capacitated Facility Location, Minimum Set Covering, Combinatorial Auctions, and Maximum Independent Set.
Randomly generated instances are solved offline using SCIP~\citep{GleixnerEtal2018OO} to collect training samples of the form $(\mathbf{G}^0, \mathbf{G}, \mathbf{X}, i^\star_{SB})$.
We leave the description of the data collection and training details of each model to the supplementary materials.

\paragraph{Evaluation.} As in \citet{gasse2019exact}, our evaluation instances are labeled as small, medium, and big based on the size of underlying MILP. Small instances have the same size as those used to generate the training datasets, and thus match the training distribution, while instances of increasing size allows us to measure the generalization ability of the trained models. Each scenario uses 20 instances, solved using 3 different seeds to account for solver variability.
We report standard metrics used in the MILP community for benchmarking \bnb{} solvers: (i) Time: 1-shifted geometric mean\footnote{for complete definition refer to Appendix A.3 in ~\citet{achterberg_thesis}} of running times in seconds, including the running times for unsolved instances, (ii) Nodes: hardware-independent 1-shifted geometric mean of \bnb{} node count of the instances solved by all branching strategies
, and (iii) Wins: number of times each branching strategy resulted in the fastest solving time, over total number of solved instances.
All branching strategies are evaluated using the open-source solver SCIP~\citep{GleixnerEtal2018OO} with a time limit of 45 minutes, and cutting planes are allowed only at the root node.

\input{inputs/model_results}
\paragraph{Baselines.} We compare our hybrid model to several non-ML baselines, including SCIP's default branching heuristic \textit{Reliability Pseudocost Branching} (RPB), the ``gold standard" heuristic \textit{Full Strong Branching} (FSB), and a very fast but ineffective heuristic \textit{Pseudocost Branching} (PB).
We also include the GNN model from \citet{gasse2019exact} run on CPU (GNN), and also several fast but less expressive models such as SVMRank from \citet{confs/aaai/KhalilBND16}, LambdaMART from  \citet{burges2010ranknet}, and ExtraTree Classifier from \citet{geurts2006extremely} as benchmarks.
For conciseness, we chose to report those last three competitor models using an optimistic aggregation scheme, by systematically choosing only the best performing method among the three (COMP).
For completeness, we also report the performance of a GNN model run on a high-end GPU, although we do not consider that method as a baseline and therefore do not include it in the Wins indicator.
We acknowledge that our comparison to general branching strategy like RPB is not completely fair, however, developing specialized versions of such strategies is a challenge in itself full of modeling choices that we foresee as future work.

\paragraph{Model selection.} To investigate the effectiveness of different architectures, we empirically compare the performance of their end-to-end variants.
Figure~\ref{fig:model_results} compares the Top-1 test accuracy of models across the four problem sets.
The performance of GNNs (blue), being the most expressive model, serves as an upper bound to the performance of hybrid models.
All of the considered hybrid models outperform MLPs (red) across all problem sets.
Additionally, we observe that FiLM (green) and CONCAT (purple) perform significantly better than other architectures.
However, there is no clear winner among them.
We also note that the cheapest hybrid models, HyperSVM and HyperSVM-FiLM, though better than MLPs, still do not perform as well as FiLM or CONCAT models.

\paragraph{Training protocols.} In Table~\ref{tabl:protocol-results}, we show the effect of different protocols discussed in section~\ref{sec:arch&train} on Top-1 accuracy of the FiLM models.
We observe that the presence of these protocols improves the model accuracy by 0.5-0.9\%, which translates to a minor yet practically useful improvement in \bnb{} performance of the solver.
Except for Combinatorial Auctions, FiLM's performance is improved by knowledge distillation, which suggests that the soft targets of pre-trained GNN yield a better generalization performance.
Lastly, we launch a hyperparameter search for auxiliary objective: ED and MHE, on top of the best performing model (Top-1 accuracy), among e2e and e2e \& KD models.
Auxiliary tasks further help the accuracy of the hybrid models, but for some problem classes it is ED that works well while for others it is MHE.
We provide the model performances on test dataset in the supplement for further reference.

\input{inputs/combined-protocol-loss-weighting}
\paragraph{Effect of loss weighting.} We empirically investigate the effect of the different loss weighting schemes discussed in Section~\ref{sec:weighting}. We train a simple MLP model on our small Combinatorial Auctions instances, and measure the resulting \bnb{} tree size on big instances. We report aggregated results in Table~\ref{tab:weights}, and provide instance-level results in the supplement.
We observe that the most commonly used exponential and linear schemes actually seem to degrade the performance of the learned strategy, as we believe those may be too aggressive at disregarding nodes early on in the tree. On the other hand, both the quadratic and sigmoidal schemes result in an improvement, thus validating the idea that depth-dependent weighting can be beneficial for learning to branch. We therefore opt for a sigmoidal loss weighting scheme in our training protocol.

\paragraph{Complete benchmark.} Finally, to evaluate the runtime performance of our hybrid approach, we replace  SCIP's default branching strategy with our best performing model, FiLM.
We observe in Table~\ref{tab:evaluation} that FiLM performs substantially better than all other CPU-based branching strategies.
The computationally expensive FSB, our ``gold standard'', becomes impractical as the size of instances grows, whereas RPB remains competitive.
While GNN retains its small number of nodes, it loses in running time performance on CPU.
Note that we found that FiLM models for Maximum Independent Set did initially overfit on small instances, such that the performance on larger instances degraded substantially.
To overcome this issue we used weight decay~\cite{ng2004feature} regularization, with a validation set of 2000 observations generated using random medium instances (not used for evaluation).
We report the performance of the regularized models in the supplement, and use the best performing model to report evaluation performance on medium and big instances.
We also show in the supplement that the cheapest computation model of HyperSVM/HyperSVM-FiLM do not provide any runtime advantages over FiLM.
We achieve up to 26\% reduction on medium instances and up to 8\% reduction on big instances in overall solver running time compared to the next-best branching strategy, including both learned and classical strategies.
\input{inputs/evaluation-compressed-format}

Finally, we note that the majority of ``big" problems in Set Covering and Maximum Independent Set are not solved by any of the branching strategies.
\input{inputs/optimality-gap}
Therefore, we provide a comparison of optimality gap of unsolved instances at time out in Table~\ref{table:opt-gap}.
It is evident that FiLM models are able to close a larger optimality gap than the other branching strategies.

In the absence of significant difference in KD and KD \& AT model's Top-1 accuracy (see Table~\ref{tabl:protocol-results}), a natural question then to ask is: what can be a practically useful choice?
To answer this question, we tested the effect of each training protocol on the final \bnb{} performance.
Although the detailed discussion is in the supplement we conclude that in the absence of significant difference in the test accuracy of models between KD and AT \& KD, a preference should be made for KD models to attain better generalization.

\paragraph{Limitations.} We would like to point out some limitations of our work.
First, given the NP-Hard nature of MILP solving, it is fairly time consuming to evaluate performance of the trained models on the instances bigger than considered for this work.
One can consider the primal-dual bound gap after a time limit as an evaluation metric for the bigger instances, but this is misaligned with the solving time objective.
Second, we have used Top-1 accuracy on the test set as a proxy for the number of nodes, but there is an inherent distribution shift because of the sequential nature of \bnb{} that leads to out-of-distribution observations.
Third, generalization to larger instances is a central problem in the design of branching strategies.
Several techniques discussed in this work form only a part of the solution to this problem.
On the Maximum Independent Set problem, we originally noticed a poor generalization capability, which we addressed by cross-validation using a small validation set.
In future work, we plan to perform an extensive study of the effect of architecture and training protocols on generalization performance.
Another experiment worth conducting would be to train on larger instances than ``small" problems used in the work, in order to get a better view of how and to which point the different models are able to generalize.
However, not only is this a time-consuming process, but there is also an upper limit on the size of the problems on which it is reasonable to conduct experiments, simply due to hardware constraints.
Finally, although we showed the efficacy of our models on a broad class of MILPs, there may be other problem classes for which our models might not result in a substantial runtime improvements.

\section{Conclusion}
As more operations research and integer programming tools start to include ML and deep learning modules, it is necessary to be mindful of the practical bottlenecks faced in that domain. To this end, we combine the expressive power of GNNs with the computational advantages of MLPs to yield novel hybrid models for learning to branch, a central component in MILP solving.
We integrate various training protocols that augment the basic MLPs and help bridge the accuracy gap with more expensive models.
This competitive accuracy translates into savings in time and nodes when used in MILP solving, as compared to both default, expert-designed branching strategies and expensive GNN models, thus obtaining the ``best of both worlds" in terms of the time-accuracy trade-off.
More broadly, our philosophy revolves around understanding the intricacies and practical constraints of MILP solving and carefully adapting deep learning techniques therein. We believe that this integrative approach is crucial to the adoption of statistical learning in exact optimization solvers.

\section*{Broader Impact}

This paper establishes a bridge between the work done in ML in the last years on learning to branch and the traditional MILP solvers used to routinely solve thousands of optimization applications in energy, telecommunications, logistics, biology, just to mention a few.

The MILP solvers are executed on CPU-only machines and the technical challenge of using GPU-based algorithmic techniques to hybridize them had been neglected thus far. Admittedly, such a challenge was not urgent when the learning to branch literature was in its early stages. That situation has changed drastically with the GNN implementation in \cite{gasse2019exact}, the first approach to show significant benefit with respect to the default version of a state-of-the-art MILP solver like SCIP. For this reason, the current paper comes at the due time for the literature in the field and addresses the challenge, for the first time, in a sophisticated, yet relatively simple way.

Thus, our work provides the first viable way for commercial and noncommercial MILP solver developers to implement and integrate a ML-based ``learning to branch" framework and for hundreds of thousands of users and practitioners to use it. In an even broader sense, the fact that we were able to approximate the performance of GPU-based models with a sophisticated integration of CPU-based techniques is  consistent with, for example, \citet{hinton2015distilling}, and widens the space of problems to which ML techniques can be successfully applied.

\begin{ack}
The authors are grateful to CIFAR and IVADO for funding and Compute Canada for computing resources.

We would further like to acknowledge the important role played by our colleagues at Mila and CERC through building a fun learning environment.
We would also like to thank Felipe Serano and Benjamin M\"uller for their technical help with SCIP and insightful discussions on branching in MILPs.
PG wants to thank Giulia Zarpellon, Didier Chételat, Antoine Prouvost, Karsten Roth, David Yu-Tung Hui, Tristan Deleu, Maksym Korablyov, and Alex Lamb for enlightening discussions on deep learning and integer programming.
\end{ack}

\small
\bibliography{references}

\end{document}


text:
\section{Introduction and main results}\label{Section1}

\subsection{Preface}\label{Section1.1}

More than thirty years ago Kardar, Parisi and Zhang \cite{KPZ} studied the time evolution of random growing interfaces and proposed the following stochastic partial differential equation (called the {\em KPZ equation}) for a height function $\mathcal{H}(T,X) \in \mathbb{R}$
\begin{equation}\label{KPZEq}
\partial_TH(X,T) = \frac{1}{2} \partial^2_X \mathcal{H}(X,T) + \frac{1}{2} [ \partial_X \mathcal{H}(X,T)]^2 +  \xi(X,T).
\end{equation}
In (\ref{KPZEq}) the letters $T$ and $X$ denote time and space, and $\xi$ is space-time Gaussian white noise, so that formally $\mathbb{E} \left[ \xi(X,T) \xi(Y,S)\right] = \delta(T-S)\delta(X-Y).$ Drawing upon the earlier work of Forster, Nelson and Stephen \cite{FNS}, KPZ predicted that for large time $T$, the height function $\mathcal{H}(X,T)$ exhibits fluctuations of order $T^{1/3}$ and has spatial correlation length of order $T^{2/3}$. The critical exponents $1/3$ and $2/3$ are believed to be universal for a large class of growth models, which is now called the KPZ universality class. For more on the KPZ universality class we refer to the surveys and books \cite{CU2, QS, HT} and the references therein.

The $3$ : $2$ : $1$ scaling of time : space : fluctuation (now known as KPZ scaling) is believed to take any model from the KPZ universality class to a universal fluctuating field, which does not depend on the particular model, but does depend on its initial data class. The conjectural space-time limiting field that attracts all models in the class is now called the {\em KPZ fixed point}, which can be described either through its Markov transition kernel \cite{MQR} or through a variational formula involving a multiparameter scaling limit known as the {\em Airy sheet} or {\em directed landscape} \cite{DOV18}. There are various natural choices for initial conditions that have been considered, see \cite[Figure 4]{CU2}, but in the present paper we will focus our attention on the {\em narrow wedge} initial data, which for discrete models is typically called {\em step} initial data. For models in the KPZ universalty class, started from step initial conditions, it is believed that under the KPZ scaling the height function converges to the {\em Airy process}. The Airy process was first introduced in \cite{Spohn} and  is a stationary continuous process, whose one point distribution is given by the celebrated Tracy-Widom distribution \cite{TWAiry}.

While the $1/3$ and $2/3$ critical exponents have been established in greater generality, demonstrating the convergence to the Airy process (in a rigorous mathematical sense) has only been done for special {\em integrable} or {\em exactly solvable} models. For the asymmetric simple exclusion process (ASEP) the one-point convergence to the Airy process, i.e. the convergence to the Tracy-Widom distribution, was first established in \cite{TWASEP}, see also \cite{TWASEP2}. Analogous asymptotic results are proved for the stochastic six-vertex model in \cite{Bor16, BCG14}, the KPZ equation in \cite{CorQ, BCF2}, the semi-discrete directed polymer in \cite{BorCor}, the log-gamma polymer in \cite{BorCorRem,KQ}, and the $q$-TASEP in \cite{Bar15,FerVet}.

Some of the aforementioned models have special cases, sometimes referred to as a {\em free fermion point} or {\em zero temperature limit} in the case of polymer models, where the multi-point (and not just single-point) convergence  to the Airy process is known. Possibly the most famous example comes from the totally asymmetric simple exclusion process (TASEP), which is obtained from the ASEP by sending the asymmetry parameter to zero. In \cite{Joh05} it was shown that when started from step initial conditions the (discrete time) TASEP, or geometric last passage percolation, converges to the Airy process in the sense of finite dimensional distributions. Analogous results exist for the polynuclear growth model \cite{Spohn, JDPNG} and random lozenge tilings \cite{FSpohn, P1}.

The common feature shared by all models for which finite-dimensional convergence to the Airy process is known is that they have the structure of determinantal processes. The determinantal structure is especially useful for proving convergence to the Airy process, since the latter is itself a determinantal process with a correlation kernel given by the {\em extended Airy kernel} \cite{Spohn}. In particular, when proving that a determinantal process converges to the Airy process in the finite dimensional sense it is sufficient to show that the kernel converges (in a sufficiently strong sense) to the extended Airy kernel, and this is indeed how convergence has been proved in the past. The problem with generalizing this approach to models like the ASEP is that once one moves away from the free fermion point, the determinantal structure is lost. The way the one-point convergence to the Airy process has been established for positive temperature integrable models (including the ASEP, stochastic six-vertex model, KPZ equation) is to utilize one of two main algebraic structures -- Macdonald processes and quantum integrable systems. In fact, there are now bridges between these structures indicating that they might be eventually joined together. Both of these structures produce moment formulas, which formally should completely characterize the distribution. However, despite the existence of multi-point exact formulas for various positive temperature integrable models \cite{BorCor, BCS, BP16, NZ} the convergence to the Airy process has proved elusive so far and only rigorously worked out in the one-point case.\\

{\bf The goal of the present paper is to prove that the height function of the stochastic six-vertex model at two points jointly converges to the Airy process, which is the first rigorous multi-point convergence result for a positive temperature integrable model in the KPZ universality class.} The starting point of our analysis is a remarkable distributional equality between the stochastic six-vertex model and the {\em ascending Hall-Littlewood process} (a special case of the Macdonald processes from \cite{BorCor}), which was established in \cite{BBW}. This identification allows us to recast the problem into the framework of Macdonald processes, where we use the method of the Macdonald difference operators from \cite{BorCor} to derive two-point observables that are suitable for asymptotic analysis. Our current framework suffers from two types of limitations: (1) we only derive formulas and perform the asymptotics for two points and (2) we can only carry out the framework for a small (but non-trivial) range of parameters. In this sense, the present paper is a proof of concept -- that one can use the method of the Macdonald difference operators to obtain multi-point convergence to the Airy process. In the future we hope to extend our framework to arbitrary parameters and number of points. In addition, we hope that the approach we develop can be extended to other integrable models in the KPZ universality class (for example by utilizing the limit of the stochastic six-vertex model to the ASEP observed first in \cite{Gwa} and proved in significant generality in \cite{Agg16}).

Despite our result being the first of its kind for the stochastic six-vertex model, there exist previous (incomplete/non-rigorous) works showing the two-point convergence of the KPZ equation started from narrow wedge initial data \cite{Dot13, Dot14, PSpohn, ISS13} and the log-gamma polymer \cite{NZ} to the Airy process. In Section \ref{Section9} we give a more detailed account of these previous works, and formally compare the techniques of the present paper with those in \cite{Dot13,PSpohn,ISS13} and \cite{NZ}. Specifically, in Section \ref{Section9.1} we give a formal explanation of how to obtain a prelimit formula using the Maconald difference operators and compare our result to the formulas obtained for the KPZ equation in \cite{Dot13,PSpohn,ISS13}. In Section \ref{Section9.2} we discuss some of the convergence issues of our formulas, and give a concise explanation of how we do the asymptotic analysis, comparing the approach to the one for the log-gamma polymer in \cite{NZ}.

Here we mention that the works addressing the KPZ equation come from the physics literature, are based on the Bethe ansatz replica technique, and at their core involve a non-rigorous moment expansion formula. As we explain in Section \ref{Section9.1}, the moment expansion formulas in \cite{Dot13,PSpohn,ISS13} can be seen as shadows of the rigorous $t$-moment formulas from the present paper; however, we are presently unable to conceptually match the Bethe ansatz framework in those papers with our difference operators approach. In \cite{NZ} the authors derive formulas for the joint Laplace transform of the partition function of the log-gamma polymer model at several locations using the geometric Robinson-Schensted-Knuth correspondence. Afterwards, using a result from \cite{BCR}, the authors are able to rewrite their formulas as a ``Fredholm determinant''-like series. These series term-wise converge to a corresponding Fredholm determinant expansion for the Airy process, and the essential ingredient missing and making the proof incomplete is an estimate on the growth of the terms in the series that would allow one to exchange the order of the sum and the limit.

The reason the authors of \cite{NZ} were unable to obtain a suitable bound on their Fredholm-like series comes from the presence of certain ``cross terms'' in the formulas that in a sense reflect the correlation of the log-gamma polymer partition functions at two locations. Part of the progress made in our paper is the ability to control similar cross terms (at least for some small range of parameters) and obtain the necessary bounds on the analogous Fredholm-like series that we derive for our model, see also Section \ref{Section9.2.2}. Nevertheless, we want to emphasize the importance of \cite{NZ} from which the present paper has greatly benefited. Indeed, as mentioned earlier for the log-gamma polymer (as is the case for all known positive temperature integrable models) the determinantal structure is lost and despite having multi-point observables it was a significant challenge to obtain any formula that would converge to the joint cdf for the Airy process. One of the many remarkable contributions of \cite{NZ} is finding a way to rewrite the joint Laplace transform for the log-gamma polymer through ingenious and highly non-trivial manipulations in a form that mimics the Fredholm determinant structure of the limit. In developing the results of the present paper, we have frequently drawn inspiration from \cite{NZ} and many of our own formulas can be seen as discrete analogues of those in \cite{NZ}.

The remainder of the introduction is structured as follows. In Section \ref{Section1.2} we define the stochastic six-vertex model on a quadrant and present our main result as Theorem \ref{thmMain}. In Section \ref{Section1.3} we give an outline of the paper and our approach. We also eagerly recommend Section \ref{Section9} to readers who are interested in a more accessible general exposition of our arguments.

\subsection{Main result}\label{Section1.2} In this section we give the definition of the homogeneous stochastic six-vertex model in a quadrant, considered in \cite{Gwa,BCG14,BP16}, and state the main result we prove about it. There are several (equivalent) ways to define the model and we follow \cite[Section 1.1.2]{Agg16B}.

A {\em six-vertex directed path ensemble} is a family of up-right directed paths drawn in the first quadrant $\mathbb{Z}^2_{\geq 1}$ of the square lattice, such that all the paths start from a left-to-right arrow entering each of the points $\{(1,m): m \geq 1 \}$ on the left boundary (no path enters from the bottom boundary) and no two paths share any horizontal or vertical edge (but common vertices are allowed); see Figure \ref{S1_1}. In particular, each vertex has six possible {\em arrow configurations}, presented in Figure \ref{S1_2}.
\begin{figure}[h]
\centering
\scalebox{0.6}{\includegraphics{S1_1.jpg}}
\caption{The left picture shows an example of a six-vertex directed path ensemble. The right picture shows an element in $P_n$ for $n = 6$. The vertices on the dashed line belong to $D_n$ and are given half of an arrow configuration if a directed path ensemble from $P_n$ is drawn.   }
\label{S1_1}
\end{figure}

\begin{figure}[h]
\centering
\scalebox{0.5}{\includegraphics{S1_2.jpg}}
\caption{The top row shows the six possible arrow configurations at a vertex. The bottom row shows the probabilities of top-right completion, given the bottom-left half of a configuration. }
\label{S1_2}
\end{figure}

The stochastic six-vertex model is a probability distribution $\mathcal{P}(b_1,b_2)$ on six-vertex directed path ensembles, which depends on two parameters $b_1, b_2$ such that
$0 < b_1,  b_2 < 1$. It is defined as the $n \rightarrow \infty$ limit of a sequence of probability measures $\mathcal{P}_n$, which are constructed as follows.

For $n \geq 1$ we consider the triangular regions $T_n = \{ (x,y) \in \mathbb{Z}_{\geq 1}^2: x + y \leq n\}$ and let $P_n$ denote the set of six-vertex directed path ensembles whose vertices are all contained in $T_n$. By convention, the set $P_1$ consists of a single empty ensemble. We construct a consistent family of probability distributions $\mathcal{P}_n$ on $P_n$ (in the sense that the restriction of a random element sampled from $\mathcal{P}_{n+1}$ to $T_{n}$ has law $\mathcal{P}_n$) by induction on $n$, starting from $\mathcal{P}_1$, which is just the delta mass at the single element in $P_1$.

For any $n \in \mathbb{N}$ we define $\mathcal{P}_{n+1}$ from $\mathcal{P}_n$ in the following Markovian way. Start by sampling a directed path ensemble $\mathcal{E}_n$ on $T_n$ according to $\mathcal{P}_n$. This gives arrow configurations (as in Figure \ref{S1_2}) to all vertices in $T_{n-1}$. In addition, each vertex in $D_n = \{ (x,y) \in \mathbb{Z}^2_{\geq 1}: x+ y = n\}$  is given ``half'' of an arrow configuration, meaning that the arrows entering the vertex from the bottom or left are specified, but not those leaving from the top or right; see the right part of Figure \ref{S1_1}.

To extend $\mathcal{E}_n$ to a path ensemble on $T_{n+1}$, we must ``complete'' the configurations, i.e. specify the top and right arrows, for the vertices on $D_n$. Any half-configuration at a vertex $(x,y)$ can be completed in at most two ways; selecting between these completions is done independently for each vertex in $D_n$ at random according to the probabilities given in the second row of Figure \ref{S1_2}. In this way we obtain a random ensemble $\mathcal{E}_{n+1}$ in $P_{n+1}$ and we denote its law by $\mathcal{P}_{n+1}$. One readily verifies that the distributions $\mathcal{P}_n$ are consistent and then we define $\mathcal{P} = \lim_{n \rightarrow \infty} \mathcal{P}_n$. \\

Given a six-vertex directed path ensemble on $\mathbb{Z}_{\geq 1}^2$, we define the {\em height function } $h(x,y)$ as the number of up-right paths, which intersect the
horizontal line through $y$ at or to the right of $x$. Our main result states that under suitable rescaling the two-point distribution of the random height function $h$ converges to the two-point distribution of the {\em Airy process}, and in order to state it we need to define the latter object.

The Airy process $A(t)$ is a continuous process on $\mathbb{R}$, which was introduced in \cite{Spohn}. We define it here by its fiite-dimensional distribution functions. Given $\xi_1, \dots, \xi_m \in \mathbb{R}$ and $\tau_1 < \cdots < \tau_m$ in $\mathbb{R}$ we define $f$ on $\{\tau_1, \dots, \tau_m\} \times \mathbb{R}$ through
$$f(\tau_j, x) = {\bf 1}_{(\xi_j, \infty)}(x) \mbox{ for $j = 1,\dots, m$}.$$
Then
\begin{equation}\label{FiniteDimAiry}
\mathbb{P} \left( A(\tau_1) \leq \xi_1, \dots, A(\tau_m) \leq \xi_m \right) = \det \left(I - fA f \right)_{L^2(\{\tau_1, \dots, \tau_m\} \times \mathbb{R})},
\end{equation}
where $A$ is the {\em extended Airy kernel}
\begin{equation}\label{AiryKernel}
A(\tau, \xi; \tau', \xi') = \begin{cases}\int_0^\infty e^{-\lambda (\tau - \tau') }Ai(\xi + \lambda) Ai(\xi' + \lambda) d\lambda \hspace{5mm} &\mbox{ if $\tau \geq \tau'$},\\ -\int_{-\infty}^0 e^{-\lambda (\tau - \tau') }Ai(\xi + \lambda) Ai(\xi' + \lambda) d\lambda \hspace{5mm} &\mbox{ if $\tau < \tau'$}.\end{cases}
\end{equation}
and $Ai(\cdot)$ is the Airy function. When $\tau = \tau'$ the extended Airy kernel reduces to the usual Airy kernel from \cite{TWAiry}. In (\ref{FiniteDimAiry}) the $L^2$ space is defined with respect to the product measure on $\{\tau_1, \dots, \tau_m\} \times \mathbb{R}$ coming from the counting measure on $\{\tau_1, \dots, \tau_m\} $ and the usual Lebesgue measure on $\mathbb{R}$. In \cite{JDPNG} it was shown that $fA f $ is a trace class operator on $L^2(\{\tau_1, \dots, \tau_m\} \times \mathbb{R})$, so that the determinant in (\ref{FiniteDimAiry}) is the usual Fredholm determinant of trace class operators, see \cite{Simon}. Numerically, the Fredholm determinant in (\ref{FiniteDimAiry}) is equal to
\begin{equation}\label{FDE}
 1+ \sum_{n = 1}^\infty \frac{(-1)^n}{n!} \sum_{i_1, \dots, i_n = 1}^m  \int_{\xi_{i_1}}^\infty \cdots \int_{\xi_{i_n}}^\infty \det \left[ A(\tau_{i_k}, x_k; \tau_{i_l}, x_l)  \right]_{k,l= 1}^n dx_n \cdots dx_1,
\end{equation}
where the latter sum converges absolutely as the Fredholm series expansion of a trace class operator.

The main result of the paper is as follows.
\begin{theorem}\label{thmMain}Let $0 < b_1 < b_2 < 1$ and put $a = \sqrt{\frac{1- b_2}{1-b_1}}$, $t = \frac{b_1}{b_2}$. There exist $a^*, t^* \in (0,1)$ such that that the following holds for any $0 < b_1 < b_2 < 1$ that satisfy $a \in (0, a^*]$, $t \in (0,t^*]$. Let $s_1, s_2 \in \mathbb{R}$ be such that $s_1 > s_2$. For $M \in \mathbb{N}$ sufficiently large so that $M + s_2 M^{2/3} \geq 1$ we define $n_1(M), n_2(M) \in \mathbb{N}$ through
\begin{equation}\label{S1ScaleN}
n_1 =  \lfloor M + s_1 M^{2/3} \rfloor \mbox{ and } n_2 = \lfloor M + s_2 M^{2/3} \rfloor.
\end{equation}
Let $h(x,y)$ be the height function of the stochastic six-vertex model distributed according to $\mathcal{P}(b_1,b_2)$. Define the random variables
$$\tilde{X}_M = \sigma_a^{-1}M^{-1/3}( - h(n_1+1 )- (f_1-1)M - f_1' [M - n_1] - (1/2)f''_1 s_1^2 M^{1/3} )$$
$$\tilde{Y}_M= \sigma_a^{-1}M^{-1/3}( - h(n_2 + 1) - (f_1-1)M - f_1' [M - n_2] - (1/2)f''_1 s_2^2 M^{1/3} ),$$
where
\begin{equation}\label{S1eqnConst}
\sigma_a = \frac{a^{1/3} \left(1 - a \right)^{1/3} }{1 + a }, \hspace{3mm} f_1=   \frac{2a}{1 + a}, \hspace{3mm} f_1'=  \frac{a}{1+a}, \hspace{3mm}  f_1'' = \frac{-a}{2  (1 - a^2)}.
\end{equation}
Then we have that for any $x_1, x_2 \in \mathbb{R}$
\begin{equation}\label{limitSSV}
\lim_{M \rightarrow \infty} \mathbb{P} \left(\tilde{X}_M \leq x_1, \tilde{Y}_M \leq x_2 \right) = \mathbb{P} \left(A(\tau_1) \leq x_1, A(\tau_2) \leq x_2 \right),
\end{equation}
where $A(\cdot)$ is the Airy process from (\ref{FiniteDimAiry}) and $\tau_1 = \frac{s_1a^{1/3}}{2(1-a)^{2/3}}$, $\tau_2 = \frac{s_2a^{1/3}}{2(1-a)^{2/3}}$.
\end{theorem}
\begin{remark}\label{S1R1} The condition $b_1 < b_2$ is necessary to obtain the Airy process limit and if $b_1 \geq b_2$ a different asymptotic behavior is expected, cf. \cite{BCG14}. In the coordinates $(b_1, b_2)$, the range $0 < b_1 < b_2 < 1$ defines an open right triangle in $\mathbb{R}^2$ with vertices at $(0,0)$, $(0,1)$ and $(1,1)$. In simple words, Theorem \ref{thmMain} states that if $(b_1, b_2)$ is close to the point $(0,1)$ then the height function of the stochastic six-vertex model evaluated at two points near the $x = y$ line, properly shifted and scaled, converges jointly to the two-point distribution of the Airy process.  The assumption that $(b_1, b_2)$ is sufficiently close to $(0,1)$ is technical,  and we believe that the theorem should hold even if we remove it. We discuss this parameter limitation later in Remarks \ref{RemS2} and \ref{RemarkRest}. The proof of the Theorem \ref{thmMain} can be found in Section \ref{Section2.2}.

We choose to evaluate the height function at $n_i +1$ in the definition of $\tilde{X}_M,\tilde{Y}_M$ (and not $n_i$) in order to obtain slightly simpler formulas later in the text. Since $|h(n_i+1) - h(n_i)| \leq 1$ we see that shifting the argument of $h$ by finite quantities does not affect the limit because of the $M^{1/3}$ scaling.
\end{remark}

\subsection{Outline}\label{Section1.3} In this section we give a brief outline of the general approach we take to prove Theorem \ref{thmMain}. The discussion below will involve certain expressions that will be properly introduced in the main text, and which should be treated as black boxes for the purposes of the outline.

The starting point of the proof of Theorem \ref{thmMain} is a remarkable distributional equality between the stochastic six-vertex model and the ascending Hall-Littlewood process, which is a measure on sequences of partitions $\Lambda = (\lambda(1), \dots, \lambda(N))$ that depends on two sets of parameters $X = (x_1, \dots, x_N)$ and $Y = (y_1, \dots, y_M)$ and a number $t\in [0,1)$. The probability of $\Lambda$ is given by
$$\mathbb{P}_{X,Y}(\Lambda) = \prod_{i = 1}^N \prod_{j = 1}^M \frac{1 - x_iy_j}{1 - t x_i y_j} \times \prod_{i = 1}^N P_{\lambda(i) / \lambda(i-1)} (x_i) \times Q_{\lambda(N)}(Y),$$
where $P_{\lambda/ \mu}$ and $Q_{\lambda}$ denote the (skew) Hall-Littlewood polynomials with parameter $t$, see \cite[Chapter 3]{Mac}. The ascending Hall-Littlewood process is defined in Section \ref{Section2.2}.

As a special case of \cite[Theorem 4.1]{BBW}, we have that if $x_i = y_j = a$ for all $i = 1, \dots, N$ and $j = 1, \dots, M$ then the following distributional equality holds
\begin{equation}\label{S1DEQ}
 \left(M - \lambda_1'(0),\dots,M - \lambda'_1(N) \right) \eqd \left(h(1,M), \dots, h(N+1,M) \right),
\end{equation}
where $\lambda'_1$ denotes the largest column of $\lambda$, $\lambda_1'(0) \equiv 0$ and $h$ is the height function of the stochastic six-vertex model distributed according to $\mathcal{P}(b_1, b_2)$ with $b_1/b_2 = t$ and $(1-b_2)/(1-b_1) = a^2$.

In view of (\ref{S1DEQ}) we see that Theorem \ref{thmMain} can be rephrased in terms of the ascending Hall-Littlewood process by simply replacing $h(n_i+1)$ with $M - \lambda_1'(n_i)$ everywhere (this is why we shifted the argument in $h$ by $1$, see Remark \ref{S1R1}). This restatement can be found as Theorem \ref{thmHLMain} in the main text. The benefit of recasting the problem in the setup of the Hall-Littlewood process is that we can apply the Macdonald difference operators to obtain joint observables for $\lambda_1'(n_1)$ and $\lambda_1'(n_2)$. We recall these operators in Section \ref{Section2.3} and eventually work with their affine shifted version that we denote by $\mathcal{D}_n$ (here $n$ means that the operator acts on the variables $x_1, \dots, x_n$).

The operator $\mathcal{D}_n$ is an eigenoperator for the Hall-Littlewood polynomials $P_\lambda(x_1, \dots, x_n)$ with eigenvalue $t^{-\lambda_1'}$ and by utilizing this fact alone one can obtain the formula
$$\mathbb{E}_{X,Y}\left[ t^{-k \lambda_1'(n_1)} \right] =  \frac{\mathcal{D}^k_{n_1} \Pi(X;Y) }{\Pi(X;Y)} , \mbox{ where }\Pi(X;Y) =\prod_{i = 1}^N \prod_{j = 1}^M \frac{1 - t x_i y_j}{1 - x_iy_j}, $$
and we have written $\mathbb{E}_{X,Y}$ for the expectation with respect to $\mathbb{P}_{X,Y}$. The expression on the right side of the above $k$-th moment formula can be written as a $k$-fold {\em nested contour integral}, where the contours become larger as $k$ increases. In order to handle this problem of growing contours, we deform all of these contours to the same one. As one deforms all of these contours, certain poles are crossed that diminish the dimension of the integral, and the nested contour integral can be rewritten as a sum over residue subspaces of integrals over the same contour. This is an instance of the {\em nested contour integral ansatz}, which was investigated in great detail in \cite{BBC}. The result of applying this ansatz method is that
\begin{equation}\label{S1Moment1}
\mathbb{E}_{X,Y}\left[ t^{-k \lambda_1'(n_1)} \right]  = \sum_{\lambda \vdash k} \frac{1}{(2\pi \iota)^{\ell(\lambda)}}\int_{\gamma^{\ell(\lambda)}} \det \left[\frac{1}{z_j t^{-\lambda_i} - z_i} \right]_{i,j = 1}^{\ell(\lambda)} \prod_{i = 1}^{\ell(\lambda)} F_{n_1}(z_i,\lambda_i; X,Y) d\vec{z} .
\end{equation}
In (\ref{S1Moment1}) the sum is over partitions of $k$, which are denoted by $\lambda$ -- these are the labels of the residue subspaces that come from the contour integral ansatz. The contour $\gamma$ is a positively oriented, zero-centered circle that contains $x_1, \dots, x_{n_1}$ and excludes $y_1^{-1}, \dots, y_M^{-1}$, $\ell(\lambda)$ is the number of parts of $\lambda$, $\iota = \sqrt{-1}$, and the function $F_{n}(z,m; X,Y) $ is given by
$$ F_{n}(z,m; X,Y) = \prod_{j =1}^M \frac{1 - zy_j}{1 - zt^{m} y_j} \cdot \prod_{j = 1}^n \frac{1 - z^{-1} t^{-m} x_j}{1 - z^{-1} x_j}.$$

Starting from (\ref{S1Moment1}) one can use the generating series of the $t$-exponential function
\begin{equation}\label{S1tExpon}
e_t(u) = \frac{1}{((1-t)u;t)_\infty} = \sum_{k = 0}^\infty \frac{u^k(1-t)^{-k}}{k_t!},
\end{equation}
where $(a;t)_\infty = \prod_{m = 0}^\infty (1 -at^m)$ is the $t$-Pochhammer symbol and $k_t! = \frac{(1 - t)(1-t^2) \cdots (1-t^k)}{(1-t)^k}$, to obtain
\begin{equation}\label{S1Laplace1}
\begin{split}
&\mathbb{E}_{X,Y}\left[ \frac{1}{(u_1t^{-\lambda_1'(n_1)};t)_{\infty}} \right] = \sum_{k = 0}^\infty \frac{u_1^k(1-t)^{-k}}{k_t!} \mathbb{E}_{X,Y}\left[ t^{-k \lambda_1'(n_1)} \right]  = \\
& \sum_{k = 0}^\infty \sum_{\lambda \vdash k} \frac{u_1^k(1-t)^{-k}}{k_t!} \frac{1}{(2\pi \iota)^{\ell(\lambda)}}\int_{\gamma^{\ell(\lambda)}} \det \left[\frac{1}{z_j t^{-\lambda_i} - z_i} \right]_{i,j = 1}^{\ell(\lambda)} \prod_{i = 1}^{\ell(\lambda)} F_{n_1}(z_i,\lambda_i; X,Y) d\vec{z} .
\end{split}
\end{equation}
In the top row of (\ref{S1Laplace1}) the expectation is a certain discrete analogue of the Laplace transform, called the $t$-Laplace transform (the name comes from the connection between the observable and the $t$-exponential function). In the second line of (\ref{S1Laplace1}) one proceeds to symmetrize the expression in $\lambda_1, \dots, \lambda_{\ell(\lambda)}$ and rewrite the sums over $\lambda$ as contour integrals, by using the formal identity
$$\sum_{n = 1}^\infty u^n g(t^n) = \frac{1}{2\pi \iota} \int_{-\iota \infty + 1/2}^{\iota \infty + 1/2} \frac{\pi}{\sin(-\pi s)} (-u)^s g(t^s)ds,$$
which essentially follows from $Res_{s = n} \frac{\pi}{\sin(-\pi s)} = (-1)^{n+1} .$ The result of this operation is that
\begin{equation}\label{S1Laplace2}
\begin{split}
&\mathbb{E}_{X,Y}\left[ \frac{1}{(u_1t^{-\lambda_1'(n_1)};t)_{\infty}} \right] = \sum_{N_1 = 0}^\infty \frac{1}{N_1!}  \frac{1}{(2\pi \iota)^{2N_1}}\int_{\gamma_1^{N_1}}\int_{\gamma_2^{N_1}} D(\vec{z}, \vec{w}) \cdot {G}_{n_1}(\vec{z}, \vec{w}; u_1, N_1) d\vec{w} d\vec{z},
\end{split}
\end{equation}
where $D(\vec{z}, \vec{w}) = \det \left[\frac{1}{z_i - w_j} \right]_{i,j = 1}^{N_1}$ is the Cauchy determinant and $\gamma_1, \gamma_2$ are two zero-centered circles with radii $r_1 > r_2$ respectively. We forgo stating what ${G}_{n_1}(\vec{z}, \vec{w}; u_1, N_1)$ is here, as it is a bit involved, but refer the interested reader to (\ref{Bu}) where the full formula is written.

We remark that in order to carry out the manipulations in (\ref{S1Laplace1}) and (\ref{S1Laplace2}) one needs to first restrict the parameters $X,Y,u_1$ to small neighborhoods of $0$ but then both sides in (\ref{S1Laplace2}) can be extended analytically to a general set of parameters. The above framework of deriving (\ref{S1Laplace2}) was carried out by the author in \cite{ED} and the formula is recalled in the main text as Lemma \ref{0thmoment}. We also mention that the right side of (\ref{S1Laplace2}) is in fact a Fredholm determinant, and equation (\ref{S1Laplace2}) was the starting point in \cite{ED}  for proving that $\lambda'_1(n_1)$ converges to the Tracy-Widom distribution.\\

In Section \ref{Section3} of the present paper we start from (\ref{S1Laplace2}) and essentially repeat the same steps above but with $\mathcal{D}_{n_2}$ instead of $\mathcal{D}_{n_1}$. There are some technical subtleties in carrying out the same framework, which will be discussed in Section \ref{Section3}, but the resulting formula has the form
 \begin{equation}\label{S1JointLaplace}
\begin{split}
&\mathbb{E}_{X,Y}\left[ \frac{1}{(u_1t^{-\lambda_1'(n_1)};t)_{\infty}} \frac{1}{(u_2t^{-\lambda_1'(n_2)};t)_{\infty}} \right] = \sum_{N_1 = 0}^\infty\sum_{N_2 = 0}^\infty \frac{1}{N_1! N_2!}  \frac{1}{(2\pi \iota)^{2N_1 + 2N_2}}\int_{\gamma_1^{N_1}}\int_{\gamma_2^{N_1}} \\
&\int_{\gamma_3^{N_2}}\int_{\gamma_4^{N_2}} D(\vec{z}, \vec{w})  {G}_{n_1}(\vec{z}, \vec{w}; u_1, N_1) \cdot D(\vec{\hat{z}}, \vec{\hat{w}})  {G}_{n_2}(\vec{\hat{z}}, \vec{\hat{w}}; u_2, N_2)  \cdot CT( \vec{z}, \vec{w}; \vec{\hat{z}}, \vec{\hat{w}}) d\vec{\hat{w}} d\vec{\hat{z}}d\vec{w} d\vec{z},
\end{split}
\end{equation}
where $\gamma_i$ are zero-centered circles with radii $r_1 > r_2 > r_3 > r_4$. As before $D(\vec{u}, \vec{v})$ stands for the Cauchy determinant and we forgo stating the exact formula for the $G$ functions, see Theorem \ref{PrelimitT} for an explicit expression of the right side in (\ref{S1JointLaplace}).

Comparing (\ref{S1JointLaplace}) with (\ref{S1Laplace2}) we see that the integrand in (\ref{S1JointLaplace}) is a product of three terms, one corresponding to $\lambda_1'(n_1)$, one corresponding to $\lambda_1'(n_2)$, but also there is now a {\em cross term} $CT$, which explicitly is given by
\begin{equation}\label{S1CT}
CT = \prod_{i = 1}^{N_1} \prod_{j = 1}^{N_2} \frac{(\hat{z}_j z_i^{-1}; t)_\infty }{(\hat{w}_j z_i^{-1} ; t)_\infty }\frac{(\hat{w}_j w_i^{-1} ; t)_\infty }{(\hat{z}_j w_i^{-1}; t)_\infty },
\end{equation}
and in a sense reflects the correlation between $\lambda_1'(n_1)$ and $\lambda_1'(n_2)$. This cross term is pointwise of order $e^{c_t N_1 N_2}$, which makes the $N_1! N_2!$ in (\ref{S1JointLaplace}) insufficient to ensure the summability of the terms. Part of the technical work behind deriving (\ref{S1JointLaplace}) in Section \ref{Section3} is showing that the Cauchy determinants in (\ref{S1JointLaplace}) provide {\em some} decay which can offset the contribution of this cross term, but we can only accomplish this if the parameters $X,Y,t$ are in a small enough neighborhood of zero, see Remark \ref{RemarkRest} for details. This is one of the sources of the parameter restriction in Theorem \ref{thmMain}. We also refer the interested reader to Section \ref{Section9.2.2} for a concise discussion on how we handle the cross term in our analysis.

Once formula (\ref{S1JointLaplace}) is established, we set all $X$ and $Y$ parameters to be equal to the same number $a$ and take the limit as $M \rightarrow \infty$. Showing that each summand on the right of (\ref{S1JointLaplace}) converges as $M\rightarrow \infty$ is an essentially straightforward application of the steepest descent argument. The exact limit statement is given as Proposition \ref{PropTermConv} and proved in Section \ref{Section5} in the main text. In order to show that the limit of the sum in (\ref{S1JointLaplace}) is equal to the sum of the limits, we require uniform in $M$ estimates on the growth (in terms of $N_1, N_2$) of the summands in (\ref{S1JointLaplace}). The order of growth we can establish is given as Proposition \ref{PropTermBound} and proved in Section \ref{Section6}. We remark that we can only find a uniform in $M$  bound on the growth of the summands if $a$ and $t$ are sufficiently close to zero -- this is the other source of the parameter restriction in Theorem \ref{thmMain}. Once Propositions \ref{PropTermConv} and Proposition \ref{PropTermBound} are proved, we know that the right side (\ref{S1JointLaplace}) converges to a certain double infinite series, which in Proposition \ref{PropTermLimit} is identified with the Fredholm determinant expansion of the two-point joint cdf of the Airy process. The weak joint convergence of $\lambda_1'(n_1), \lambda_1'(n_2)$ to the Airy process is an easy consequence of Propositions \ref{PropTermConv}, \ref{PropTermBound} and \ref{PropTermLimit}, and the argument is the content of Section \ref{Section4.2}.

\subsection*{Acknowledgments}\label{Section1.5} The author would like to thank Alexei Borodin, Guillaume Barraquand and Ivan Corwin for useful comments on earlier drafts of this paper as well as Amol Aggarwal for stimulating conversations about computing $L^2$ norms of Cauchy determinants. The author is partially supported by the Minerva Foundation Fellowship.

config: arxiv


TITLE: if $a_{n+1}=xa_{n}+ya_{n-1}$ then have this results :$a_{m+n}=Xa_{m}a_{n+1}+Ya_{m-1}a_{n}$？
QUESTION [2 upvotes]: it is well  known:
if $F_{1}=1,F_{2}=1$,and $F_{n+1}=F_{n}+F_{n-1}$.then we have
$$F_{m+n}=F_{m}F_{n+1}+F_{m-1}F_{n}$$
prove 1
and I have know this :if $a_{0}=0,a_{1}=1$ and $x,y$ be give postive integers,and $$a_{n+1}=xa_{n}+ya_{n-1}$$
then we have
$$a_{m+n}=a_{m}a_{n+1}+ya_{m-1}a_{n}$$
prove 2
My question is
if the sequence such $a,b,x,y$ be give postive integers,and such that  $a_{1}=a,a_{2}=b,a,b\in N^{+}$,and $$a_{n+1}=xa_{n}+ya_{n-1}$$,then is also have this  Similar results
$$a_{m+n}=Xa_{m}a_{n+1}+Ya_{m-1}a_{n}?$$
if there exist,Find the $X,Y$?Thanks

REPLY [2 votes]: So, the answer is indeed, yes, such a relation exists (with some caveats), and we can prove it via induction, just as in the 1st answer you have linked in the question.
MAJOR EDITS: My original answer was incorrect, as pointed out by @mathlove. Here is a revised answer.
First, let us prove the induction step. For it, assuming the relation is true upto $(m,n)$, we have to prove it for $(m+1,n)$ and $(m,n+1)$, just as mentioned in the $1^{st}$ answer.
For $(m+1,n)$, we have
\begin{align}
a_{m+1+n} &= xa_{m+n}+ya_{(m-1)+n} \\
&= x(Xa_ma_{n+1}+Ya_{m-1}a_n) + y(Xa_{m-1}a_{n+1}+Ya_{m-2}a_n) \\
&= Xa_{n+1}(xa_m+ya_{m-1}) + Ya_n(xa_{m-1}+ya_{m-2}) \\
&= Xa_{m+1}a_{n+1} + Ya_ma_n
\end{align}
Similarly we can do so for $(m,n+1)$. Thus, the inductive step is proved.
For the base case, as the induction relation involves $m-1$ and $n$, the lowest values that these can have is $1$ (As $a_0$ is not defined). So, the base case for the induction must be for $(m,n) = (2,1)$.
For the base case to hold, the following relation must hold -
\begin{align}
a_{2+1} &= Xa_2a_{1+1} + Ya_{2-1}a_1 \\
\implies a_3 &= Xb^2 + Ya^2
\end{align}
From the given recursion relation, we get
\begin{align}
a_3 &= xa_2 + ya_1 \\
&= xb + ya
\end{align}
Hence, the given relation holds iff $X,Y$ satisfy
$$Xb^2 + Ya^2 = xb + ya$$
Addendum: I added this as this was part of my original answer, and I wanted to show where I originally went wrong.
Now, let us determine what restrictions $X$ and $Y$ must have, if such a result was to exist.
Let $m=2$. Then, the induction gives
\begin{align}
a_{n+2} &= Xa_2a_{n+1}+Ya_1a_n \\
&= (Xb) a_{n+1} + (Ya) a_n
\end{align}
Now, we have already been given that $a_{n+2}=xa_{n+1}+ya_{n}$. Thus, we get
\begin{align}
(Xb) a_{n+1} + (Ya) a_n &= xa_{n+1}+ya_{n}\\
\implies (Xb-x)a_{n+1} &= (y-Ya)a_n
\end{align}
Now, if either $X = \frac{x}{b}$ or $Y = \frac{y}{a}$, then this relation shows that the other must be zero as well. Hence, we have either $(X,Y) = \left(\frac{x}{b}, \frac{y}{a}\right)$, or
\begin{align}
a_{n+1} &= \frac{y-Ya}{Xb-x}a_n \\
\implies a_n &= \left(\frac{y-Ya}{Xb-x}\right)^{n-1} a_1
\end{align}
Note: We can easily see that in both the cases, $Xb^2 + Ya^2 = xb + ya$ holds (as it should).


TITLE: Find all 2 by 2 matrices that are orthogonal, real, and also symmetric.
QUESTION [0 upvotes]: The question is taken from PSet 6.4, Q27 of Strang's Introduction to Linear Algebra, 5th edition.

Find all 2 by 2 matrices that are orthogonal and also symmetric. Which two numbers can be eigenvalues of those two matrices?

I have assumed in what follows that he also meant to say real matrices, as I can find examples not included in the answer that are complex and also satisfy the conditions.
The answer states that orthogonality and symmetry mean that $|\lambda| = 1$, (Which is only true if the matrix is orthogonal and real, I think?) and that $\lambda$ must be real (Since the matrix is symmetric and real), so $\lambda = \pm 1$.
He then goes on to say that the matrix must be either $\pm I$, or
$
S = Q\Lambda Q^T =
\begin{bmatrix}
\cos{\theta} & -\sin{\theta}\\
\sin{\theta} & \cos{\theta}
\end{bmatrix}
\begin{bmatrix}
1 & 0\\
0 & -1
\end{bmatrix}
\begin{bmatrix}
\cos{\theta} & \sin{\theta}\\
-\sin{\theta} & \cos{\theta}
\end{bmatrix}
= \begin{bmatrix}
\cos{2\theta} & \sin{2\theta}\\
\sin{2\theta} & -\cos{2\theta}
\end{bmatrix}
$.
This seems to come from nowhere to me. It's clear that a matrix with such properties must have the form $S = Q\Lambda Q^T$, but the rotation matrix as the eigenvector matrix doesn't seem to have a motivation to me other than the fact that it just works.
I'm also not at all sure what the justification is for saying this is the only form such a matrix can take. I can't see anything wrong with it, so I'm sure this describes at least some of the required matrices, but why is it the case that it describes them all?

REPLY [0 votes]: It is perhaps worth observing that this question generalizes.  First, we can let $n>2$.  Second, we can consider complex matrices (instead of real ones) by changing symmetric to be Hermetian and changing orthogonal to be unitary.  If you have not seen these terms before, if we let $A^*=\overline(A^T)$ be the complex conjugate of the transpose of $A$, then we change $A=A^T$ to $A=A^*$ and $A^{-1}=A^T$ to $A^{-1}=A^*$
Being Hermetian forces all the eigenvalues to be real.  Being unitary forces all the eigenvalues to have norm $1$.  Together, they forces all the eigenvalues to be $\pm 1$.
Then, the spectral theorem says not only that your matrix is diagonalizable, but actually that it can be diagonalized by a a unitary matrix.  There isn't quite as nice a description of things as rotation matrices, but essentially, if we pick $n$ (unitarily) orthogonal directions and $n$ choices of sign, that will give us a matrix.  However, a lot of these choices will yield the same matrix.  For example, if all our signs are positive, then whatever directions we take, we will get the identity matrix, and if all our signs are negative, then whatever directions we take, we will get $-I$.
It turns out that the matrix will be completely determined by picking one subspace of $\mathbb C^n$ to be the eigenspace of the eigenvalue $1$.  Its orthogonal complement will be the eigenspace of the eigenvalue $-1$, and with this information, there will be a unique Hermetian matrix with those eigenspaces.
With a $2\times 2$ real matrix, you can only end up with $\pm I$ or a reflection about a line.  more generally, you end up with a "reflection about a subspace".
Let $V$ be the eigenspace of $1$, and $A_V$ the corresponding matrix.  Let $P_V=(I+A_V)/2$.  Then $P_V$ will be the orthogonal projection onto $V$, and $A_V=2PV-I$.  Of course, I don't know if you find orthogonal projections easier to think about, but at least the option is there for you.