\section{Supplementary Material}

The following sections provide additional details on our experiments, our implementations, and the tasksets generated with our method.

\subsection{Experimental Details}

\subsection{Controlling Information Overlap}
For results reported in Figure~\ref{fMIAcc} we train the Multiclass network for 200 epochs on the train classes, and measure test accuracy over 600 test tasks with the weights obtained at the best validation epoch.
For optimization, we use Adam with a tuned learning rate and other hyper-parameters set to the default values in PyTorch.
We do not use weight decay regluarization nor decay the learning rate, and set the mini-batch size to 64.

For the PCA projections of Figure~\ref{fPCAClustering}, we use the PCA implementation in scikit-learn with default hyper-parameters and simply set \texttt{ncomponents=2}.
When reporting $\vert \vert \mu_\text{train} - \mu_\text{test} \vert\vert$, we measure the Euclidean distance in the original embedding space not the projected one.

Since we do not need to optimize their parameters, we use Gaussians to measure correlation between divergence and accuracies in Table~\ref{tMeasureOverlap}.
We estimate the means and covariance from the class embeddings, and use a small diagonal damping factor of 0.001 to ensure numerical stability.

\subsection{Comparing Semantic vs Embedding Clusters}
For Table~\ref{tAccuracyComparison}, the label \emph{Easy} corresponds to $D = 0.04$, \emph{Hard} to $D = 0.96$.
For \emph{Random}, we randomly partition classes among train, validation, and test sets which is equivalent to $D = 0$.
The CNN4 and ResNet12 implementeations are provided with our code release.

When comparing clustering trees in Table~\ref{tHopDistances}, we use the agglomerative clustering implementation found in scikit-learn, and set the linkage argument to \texttt{'ward'}.
We parse the structures in the following 2 files to build the WordNet graph over ImageNet-1k:
\begin{itemize}
    \small
    \item \url{http://www.image-net.org/api/xml/ReleaseStatus.xml} and
    \item \url{http://www.image-net.org/api/xml/structure_released.xml},
\end{itemize}
and prune the graph to discard leaves irrelevant to ImageNet-1k synsets.

\subsection{Is a Good Embedding Really Enough?}
For each method, the slopes in Table~\ref{tMethodSlope} are computed using the test accuracies when varying $R$ from 0.04 to 0.96.
In Table~\ref{tHardPartitions}, we report the test accuracies when $D = 0.96$.

\subsection{Code Release}
With this submission, we include implementation for \ourmethod, as well as parallelized implementations for popular few-shot learning algorithms.
(MAML, ANIL, and ProtoNet)

The few-shot learning implementations proceed one task at a time.
Before processing the task, the model is replicated across all GPUs; then, the data of the task is split across GPUs and processed by each of the replicas independently.
(This is parallelization is made simple thanks to PyTorch's \texttt{DataParallel} wrapper.)
Finally, the output of each replica is gathered on the first GPU to compute the cross-entropy loss.
Those implementations are available in \texttt{maml.py}, \texttt{anil.py}, and \texttt{protonet.py}.

The \ourmethod implementation is provided in \texttt{atg.py}.
It (and our other implementations) can be ran as-is, using the python command:
\begin{verbatim}
    python cvpr_code/atg.py
\end{verbatim}
For more details on our implementations, please refer to the included \texttt{README.md} file.
As mentioned in the main text, our implementations are available at

\noindent
\url{\codelink}.

\subsubsection{Software Acknowledgement}
This manuscript is made possible thanks to the following open-source software: Python~\cite{python}, Numpy~\cite{numpy}, PyTorch~\cite{NEURIPS2019_9015}, torchvision~\cite{torchvision}, scikit-learn~\cite{scikit-learn}, learn2learn~\cite{Arnold2020-ss}, and matplotlib~\cite{Hunter:2007}.

\subsection{Taskset Partitions}
We describe the tasksets generated by our method for different degrees of difficulty, obtained for different values of $R = 0.04, 0.32, 0.64, 0.96$.
The following tables provide the train, validation, and test class partitions for
\begin{itemize}
    \item Table~\ref{tCifar100Partitions} for CIFAR100 as implemented in torchvision,
    \item Table~\ref{tMIPartitions} for mini-ImageNet as implemented in learn2learn,
    \item Table~\ref{tEmnistPartitions} for EMNIST as implemented in torchvision,
    \item Table~\ref{tLFW10Partitions} for LFW10 as implemented in scikit-learn, and
    \item Tables~\ref{tTIPartitions} and \ref{tTIPartitions2} for tiered-ImageNet as implemented in learn2learn.
\end{itemize}
All those partitions are obtained with \ourmethod, with the following implementation details.
First, we normalize the class embeddings to unit norm; then, optimize $\mu_\text{train}$ and $\mu_\text{test}$ to minimize the objective $J$.
To that end, we use SGD with a $0.1$ learning rate and momentum set to $0.9$ for 7,000 iterations.

\input{results/table_cifar100_partitions}
\input{results/table_mi_partitions}
\input{results/table_emnist_partitions}
\input{results/table_lfw10_partitions}
\input{results/table_ti_partitions}

\section{Background and Notation}
\label{sBackground}

We now present the few-shot learning setting, introduce notation, and review the few-shot algorithms of interest to the remainder of this paper.

\paragraph{Datasets, Tasksets, and Tasks}
We designate by \emph{dataset} a set of input-class pairs $(\mathbf{x}, y)$, where the dependent variable $y$ takes values from a finite set of classes $y_1, \dots, y_M$.
This dataset induces a classification \emph{task}, which consists of finding a predictive function mapping $\mathbf{x}$ to $y$.
One common approach to construct a set of such tasks, \emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot a \emph{taskset}, is to first sample a subset of $N$ classes, and then sample $K$ input-class pairs for each of those classes.
This setting is usually referred to as $N$-ways $K$-shots.

As we often wish to evaluate the transferability of few-shot models from one taskset to another, we first partition the $M$ classes of the dataset into train and test splits and construct tasksets from those splits.
As mentioned in Section~\ref{sIntroduction}, partitioning can either be random, in which case train and test classes come from the same distribution, or it can leverage additionaal human knowledge.
This additional information typically defines a hierarchy over the classes (\emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot semantic class relationships), which can be used for partitioning.
For example, the 100 classes of CIFAR100 are grouped into 20 superclasses, and FC100 uses 12 of those superclasses (60 classes) for training and 4 (20 classes) for validation and testing, each.
Once tasksets are defined, the goal of a few-shot learning algorithm is to find a model of the train tasks that generalizes to the test tasks.

To achieve this goal, we associate to each task $\tau$ a loss $\mathcal{L}_\tau$ and a parameterized model $p_\tau(y \mid \mathbf{x})$.
Without loss of generality, we write this model as the composition of a linear layer $w$ and a feature extractor $\phi$ (\emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot a neural network).
Then, finding a maximum likelihood solution for $\phi$ and $w$ boils down to minimizing the average task loss $\Exp[\tau]{\mathcal{L}_\tau(\phi, w)}$ over the training tasks.

\paragraph{Algorithms}
Gradient-based methods, and MAML in particular, compute the likelihood after adapting $\phi$ and $w$ for one (or a few) gradient step.
For the $i$th class, this likelihood takes the form:
\begin{align*}
    \tag{MAML}
    p_\tau (y=y_i \mid \mathbf{x}) &= \text{softmax}( \phi^\prime (\mathbf{x})^\top w_i^\prime) \\
    \text{s.t.} \quad w^\prime &= w - \alpha \nabla_{w} \mathcal{L}_\tau \\
                   \phi^\prime &= \phi - \alpha \nabla_\phi \mathcal{L}_\tau
\end{align*}
where $\alpha$ denotes the adaptation learning rate, and $w_i^\prime$ is the $i$th column of the adapted linear layer.
In the adaptation step, the gradient of the loss is computed over the set of few-shot examples; differentiating it requires Hessian-gradient products which makes MAML an expensive method for large feature extractor $\phi$.
For this reason, MAML is rarely used with architectures larger than the standard 4-layer CNN.
To scale up MAML to the 12-layer ResNet of our experiments, we implement a data-parallel version of MAML where different GPUs are responsible for different subsets of the few-shot inputs.
Refer to the Supplementary Material for more details.

To alleviate this computational burden, the authors of ANIL suggest that adapting $w$ is sufficient to claim many of MAML's benefits:
\begin{align*}
    \tag{ANIL}
    p_\tau (y=y_i \mid \mathbf{x}) &= \text{softmax}( \phi (\mathbf{x})^\top w_i^\prime) \\
    \text{s.t.} \quad w^\prime &= w - \alpha \nabla_{w} \mathcal{L}_\tau.
\end{align*}
Contrasting with MAML, the difference lies in the feature extractor $\phi$ being shared -- \emph{not} adapted -- across tasks.
In that sense ANIL resembles metric-based based methods, despite its original motivation of approximating MAML.
Since $\phi$ is responsible for the bulk of the computation, ANIL becomes increasingly more efficient as the feature extractor grows in size.\footnote{In our experiments, we saw speed-ups as large as 9.25x over MAML.}

A representative metric-based algorithm, ProtoNet also shares $\phi$ across tasks and replaces the final linear classifier by a nearest neighbour one.
It measures this nearest neighbour distance between the query embedding $\phi(\mathbf{x})$ and the $i$th class embedding $\phi_i$:
\begin{align*}
    \tag{ProtoNet}
    p_\tau (y=y_i \mid \mathbf{x}) &= \text{softmax}\left( - d(\phi (\mathbf{x}), \phi_i)\right)
\end{align*}
where $\phi_i = \Exp[\mathbf{x} \sim p_\tau(\cdot \mid y_i)]{\phi (\mathbf{x})}$ is the average embedding of the few samples with class $y_i$, and $d(\cdot, \cdot)$ is a distance function -- common choices include the Euclidean norm or the (negative) cosine similarity.

Since the ProtoNet classifier is nonparametric, it can be used with any feature extractor $\phi$.
Thus, a simple yet effective baseline, \emph{Multiclass}, consists of collapsing all train tasks into a single large tasks and learning a feature extractor via standard one-vs-rest multiclass classification.
At test-time, we can readily use this learned feature extractor with the ProtoNet classifier to solve unseen tasks.

\section{Conclusion}
\label{sDiscussion}
This manuscript focuses on the evaluation and analysis of few-shot classification methods.
Accordingly, we propose \ourmethod to generate tasksets of desired transfer difficulty, and with no requirement for additional human information about classes or their relationships.
After empiricaly validating \ourmethod, we generate tasksets to study the two main families of few-shot classification algorithms: gradient-based and metric-based.
As opposed to recent work suggesting that a good feature extractor might be enough for few-shot classification~\cite{Tian2020-gv, Raghu2019-ff}, we find that gradient-based methods outperform metric-based ones, especially when transfer is challenging.

Although seemingly contradicting, we believe the two hypotheses are compatible: when applied out of domain the metric implicitly learned by ProtoNet and Multiclass will likely require adaptation, \emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot to adjust for domain shift or to discover new features.
For similar reasons and as highlighted in our experiments, ANIL's approximation of MAML will break since it does not adapt its feature extractor to avoid expensive second-order derivatives.

On the other hand, those methods don't have to pay the price of adaptation when knowledge transfer is sufficient to reach good performance, and can thus be scald to much larger datasets.
In fact, when test tasks are similar to training, MAML's adaptation by gradient descent can even lead to overfitting especially when working with limited labelled data.

We hope our contribution can help researchers analyse few-shot learning methods and answer some of the above questions.
Future research directions for \ourmethod include extensions to regression and multi-label classification tasks.

\section{Experiments}
\label{sExperiments}

Our experiments focus on three aspects.
First, we empirically validate our proposed method to control the information overlap between train and test classes.
To provide an intuitive picture of our method, we also perform an ablation study on the different measures of information and visualize low-dimensional projections of the resulting partitions.
Second, we compare the partitions obtained from our method to the ones obtained using semantic class relationships.
Our results suggest that both types of partitions are significantly different; moreover, they confirm that our method is capable of generating partitions that are more challenging than the ones resting on semantics.
Finally, we leverage our partitions to compare gradient-based and metric-based classification methods.
Those experiments indicate that metric-based methods -- which were recently thought to be sufficient for few-shot learning -- tend to underperform in challenging transfer settings when compared to their gradient-based counterpart.

For additional experimental details, including the exact partitions for all tasksets, please see the Supplementary Material.

\subsection{Setup}
Our experiments build on architectures and datasets widely used in the computer vision and few-shot learning literature.
We denote by CNN4 the 4-layer CNN with 64 hidden described in~\cite{Snell2017-kv}, which we use for few-shot learning experiments on FC100, CIFAR-FS, EMNIST, and LFW10.
(The latter two datasets are described below.)
On mini-ImageNet and tiered-ImageNet, we use the 12-layer residual network (ResNet12) described in~\cite{Mishra2017-vm}, and add the DropBlock layers proposed in~\cite{Lee2019-cl} for regularization.
(But unlike Lee~\emph{et al.} we keep the final average pooling layers.)
For experiments involving Imagenet-1k -- and including Birds, Planes, Flowers, Fungi --  we use 3 architectures: a 121-layer DenseNet (DenseNet121;~\cite{Huang2016-qc}), a 18-layer residual network (ResNet18;~\cite{He2015-qd}), and a GoogLeNet (GoogLeNet;~\cite{Szegedy2014-rn}).
For all models, we define the feature extractor to be the architecture up to the last fully-connected layer.

All tasksets built with \ourmethod follow an identical recipe.
We pretrain $\phi$ on the same data used to compute the class embeddings $\phi_i$.
We minimize $J$ by gradient descent with $\lambda = 1$, and, unless specified otherwise, use the symmetrized version of the KL divergence for $D$.

With this manuscript, we also contribute tasksets for two datasets previously not used in the few-shot learning setting: EMNIST and LFW10.

\paragraph{EMNIST}
The Extended MNIST dataset (EMNIST;~\cite{Cohen2017-vw}) is a variant of the MNIST dataset consisting of 814,255 grayscale images of 62 handwritten characters (digits, lowercase and uppercase alphabet).
Each image contains a single character scaled to 28x28 pixels.
Using \ourmethod, we partition its classes in 37 characters for training, 12 for validation, and 13 for testing, filling a niche in the few-shot dataset landscape: EMNIST is lightweight (roughly 610Mb) with few training classes, but plenty of data ($>$10k) available per class.

\paragraph{LFW10}
We bootstrap LFW10 from the \emph{Labelled Faces in the Wild} dataset~\cite{Huang2014-hw}, which contains 13,233 pictures of 5,749 famous personalities.
Of those 5k personalities, we select all 158 that have at least 10 images in the dataset -- fewer would constrain the few-shot learning setup -- and partition them in tasksets of 94 train, 32 validation, and 32 test classes.
Each image is rescaled to 62x47 pixels with RGB colors.
LFW10 is an example of dataset for which collecting class relationships is difficult if not impossible, as it requires putting a semantic hierarchy over each individual in the dataset.

\subsection{Controlling Information Overlap}\label{sControllingInformationOverlap}

\input{results/figure_mi_acc}

We first test whether our method is effective in producing tasksets of varying difficulty.
To that end, we use \ourmethod to generate tasksets of increasing difficulty on the classes of 5 datasets: CIFAR100, mini-ImageNet, tiered-ImageNet, LFW10, and EMNIST.
First, we train a convolutional network over the entire set of classes, using the standard cross-entropy minimization.
We use a 121-layer DenseNet for CIFAR100, mini- and tiered-ImageNet dataset, and the CNN4 for LFW10 and EMNIST.
We then remove the last fully-connected layer of the network, and use this feature extractor to compute the mean embedding of each class.
Finally, we create train, validation, and test class partitions for various target divergence between train and test classes.

We measure transfer difficulty with the 5-ways 5-shots classification accuracy of a model trained with Multiclass on the train taskset only.
We use the ResNet12 for ImageNet-based tasks, and the CNN4 for the others.
Figure~\ref{fMIAcc} reports train, validation, and test accuracies as a function of the divergence for all datasets.
Across all datasets test accuracy decreases as the divergence between train and test class distributions increases.
Thus, we conclude that our method is effective in finding partitions of desired transfer difficulty.
We also observe that validation and test accuracies similarly challenge the transfer ability of the Multiclass model, which is expected as both sets of classes come from the same distribution.

\input{results/figure_pca_clustering}

\paragraph{PCA Visualizations}
We perform two more ablative studies to provide further understanding of \ourmethod.
In Figure~\ref{fPCAClustering}, we plot the 2-dimensional PCA projection of the CIFAR100 mean class embeddings, together with their assignments and the assigned centroid.
We observe that centroids separate further as the divergence increases from 0.04 to 1.28 -- the greater the divergence between $p_\text{train}$ and $p_\text{test}$, the further $\mu_\text{train}$ and $\mu_\text{test}$.

\paragraph{Comparing Divergences}
Although intuitive, this description is not sufficient to explain why \ourmethod works because it lacks a measure of how difficult it is to discriminate between the classes of a task.
To show this, we compare different methods for measuring the distance of train and test class distributions between ImageNet-1k and 4 downstream datasets.
(Birds, Planes, Fungi, and Flowers)
We take a 121-layer DenseNet pretrained on ImageNet-1k, and compute the mean class embedding for all 1,000 classes as well as the ImageNet centroid.
Then, for each of the downstream dataset, we sample 100 tasks, each consisting of 5 classes and 5 samples.
Finally, we compute the distance between the train distribution (defined by the ImageNet centroid) and the task distribution (defined by the centroid of the task) over all 1,005 classes.

\input{results/table_measure_overlap}

Table~\ref{tMeasureOverlap} reports the Pearson’s correlation coefficient between those distance measurements and the task accuracies for all datasets.
The symmetrized KL, which we use in our remaining experiments, performs best while naively measuring the Euclidean distance between centroids (Euclidan) performs worse.
As a final remark, let us observe that correlations have fairly low magnitudes, indicating that our method is ill-suited to accurately measure individual task difficulty.
Assessing task similarity and difficulty is an open research question; we refer the reader to related litterature~\cite{Achille2019-hu, Fifty2020-my, Dhillon2019-mg, Wu2019-kw} for more details.

\subsection{Comparing Semantic vs Embedding Clusters}
We continue the study of our method and zero-in on the difference between class embeddings and semantic to encode class similarities.

\paragraph{Comparison to WordNet}
The first question we ask is whether class embeddings capture information similar to semantics.
To that end, we compare hierarchies created by clustering class embeddings to WordNet, the hierarchy over ImageNet-1k classes induced by class semantics.

We construct three sets of ImageNet-1k embeddings with three different pretrained architectures.
For each set, we obtain a tree of their classes via Ward clustering.
To compare those trees against the WordNet graph of classes, we define the \emph{hop} distance -- the average difference of distances between two classes $a, b$:
\begin{equation*}
    \sum_{a, b} \vert d_\text{Clustering}(a, b) - d_\text{WordNet}(a, b) \vert,
\end{equation*}
where $d_\text{Clustering}, d_\text{WordNet}$ are the normalized minimum number of nodes separating $a$ to $b$ in the clustering tree and WordNet graph, respectively.
We use such a cumbersome metric because WordNet contains cycles as well as several vertices of degree 2; intuitively, this metric can be understood as the average difference in the number ``hops’’ necessary to reach $b$ starting from $a$ between the two graphs.

\input{results/table_hop_distances}

Table~\ref{tHopDistances} reports the hop distance between WordNet and clustering trees, when the class embeddings are computed using DenseNet121, GoogLeNet (embedding sizes 104) and ResNet18. (embedding size 512)
Regardless of embedding size or network architecture, the clustering structures differ from the WordNet structure by a factor of 4, indicating that class embeddings encode attributes significantly different from class semantics.

Anecdotally, we manually inspect the trees generated from embeddings in an attempt to shed some light on what those poperties might be.
We find that the mean class embeddings tended to encode visual properties.
For example, the \emph{traffic light} and \emph{theatre} classes are clustered close-by (3 hops) due to both containing many red pixels.
On the other hand, the Orangutan and Macaque classes (fairly close semantic-wise) are separated by 9 hops — likely due to the difference in color and texture of both animal’s furs.

\paragraph{Comparison across Methods}
Having established that \ourmethod can generate tasksets different from those built on semantics, we answer the next natural question: \emph{How do those tasksets compare to existing tasksets?}
Accordingly, we select the easy ($D=0.04$) and hard ($D=0.96$) tasksets generated with \ourmethod and train networks with 4 popular few-shot methods (MAML, ANIL, ProtoNet, Multiclass) on the 5-ways 1-shot and 5-ways 5-shots setting on the tasksets of CIFAR100, mini-ImageNet, and tiered-ImageNet.
To ensure fair comparison, we do not augment the data and use the exact same architecture for all methods within a dataset.
(CNN4 for CIFAR100, ResNet12 for mini/tiered-ImageNet.%
\footnote{Training a ResNet12 with MAML on tiered-ImageNet was \emph{still} too computationally expensive, despite our data-parallel implementation.})
Each (method, taskset) pair is tuned independently, using a logarithmic-spaced grid-search.
For MAML and ANIL, we measure accuracy after 5 adaptation steps.

\input{results/table_accuracy_comparison}

Table~\ref{tAccuracyComparison} reports the test accuracy obtained at the best validation iteration.
In almost all scenarios, \ourmethod is able to generate tasksets that are as — or even more — challenging as the tasksets constructed from semantic relationships.
This makes \ourmethod a compelling solution to taskset generation in settings where additional inter-class information is hard or even impossible to obtain.
Conversely, we see that the easier tasksets generate by our method (the ones with 0.04 divergence) are significantly easier, sometimes even approaching purely random assignments. (\emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot with a divergence of 0)

\subsection{Is a Good Embedding Really Enough?}

\input{results/table_method_slope}

A closer inspection of the results in Table~\ref{tAccuracyComparison} hints at an unexpected trend: \emph{MAML becomes increasingly more competitive as the transfer from train to test classes become more challenging.}
To verify this hypothesis, we train all 4 few-shot algorithms on all tasksets of 4 datasets highlighted in Section~\ref{sControllingInformationOverlap}; then, we compute the slope of the regression line between divergence and test accuracies.
Table~\ref{tMethodSlope} reports those slopes, and confirms our hypothesis: gradient-based methods degrade slower than metric-based ones.

In turn, those experiments suggest a follow-up hypothesis: \emph{when train-test transfer is challenging enough, methods that adapt their embedding function should outperform those that do not}.
We verify this hypothesis in Table~\ref{tHardPartitions} where we compare accuracies of each algorithms on the most challenging taskset of each dataset.
We also include a \emph{Finetune} entry, which corresponds to Multiclass with the embedding function updated by 5 gradient steps on test tasks.
Confirming our hypothesis, MAML and Finetune dominate on all datasets.

\section{Introduction}
\label{sIntroduction}

Few-shot learning, the ability to learn from limited supervision, is essential for the real-world deployment of adaptive machines.
Although proposed more than 20 years ago~\cite{Miller2000-hy, Fei-Fei2006-od}, this field has recently been the focus of vast research efforts and a plethora of methods were proposed to tackle many of its challenges, including knowledge transfer and adaptation.

However, the fundamental problem of \emph{evaluating} those methods remains largely unaddressed.
Although many standardized benchmarks exist, they follow one of two recipes to generate classification tasks -- they either partition classes at random (\emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot~\cite{Bertinetto2018-gn, Vinyals2016-cm}) or leverage class semantic relationships (\emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot~\cite{Oreshkin2018-yi, Ren2018-hk, Triantafillou2019-bj}).
The former implicitly assumes that train and test tasks come from the same distribution, leading to overly optimistic evaluation.
The latter, although more realistic, requires additional human knowledge which can be expensive to gather, when available at all.
This status quo is unsatisfactory because different applications call for different benchmarking schemes: a model that performs best when train and test tasks are similar does not necessarily achieve top accuracy when the two tasksets significantly differ.
In other words, the quality of a few-shot learning algorithm depends on both train and test tasks, and \emph{their relative similarity}.

\input{results/table_hard_partitions}

Without more fine-grained benchmarks, we might miss important properties of our algorithms and hinder their deployment to real-world scenarios.
For example, recent work suggests that simply learning a good feature extractor might be \emph{all we need} for few-shot classification~\cite{Tian2020-gv, Raghu2019-ff}.
Notably, they match and often surpass the performance of gradient-based methods by sharing a feature extractor across all tasks and only adapting a final classification layer.
But, it stands to reason that in the extreme case where train tasks contain no information relevant to the test ones (\emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot transfer is impossible), those methods will underperform those that are allowed to adapt their feature extractor.
Table~\ref{tHardPartitions} diplays results for such an instances where gradient-based methods (MAML~\cite{Finn2017-gw}) dominate on tasksets carefully designed to challenge transfer.

\input{results/figure_pipeline}

In this manuscript we propose \emph{Automatic Taskset Generation} (\ourmethod), a method to automatically generate tasksets from existing datasets, with fine-grained control over the transfer difficulty between train and test tasks.
(\emph{c.f}\onedot} \def\Cf{\emph{C.f}\onedot Figure~\ref{fPipeline})
Our method can be understood as a penalized clustering objective that enforces a desired similarity between train and test tasks.
Importantly, it does not require additional human knowledge and is thus amenable to settings where this information is not available.
We use \ourmethod to study and evaluate the two main families of few-shot classification algorithms: gradient-based and metric-based methods.
Our results on 5 tasksets, including two new ones, show that gradient-based method become particularly compelling when transfer is most challenging.

\paragraph{Contributions}
We make the following contributions:
\begin{itemize}
    \item \ourmethod, a method to automatically generate tasksets that does not require additional human knowledge.
    \item Extensive validation and study of our method, showing it can effectively control the degree of transfer difficulty between train and testk tasks.
    \item An empirical analysis of popular few-shot learning methods, suggesting that gradient-based method outperform metric-based ones in the most challenging transfer regimes.
\end{itemize}
Our implementations and tasksets are described in the Supplementary Material, and are publicly available at:

\noindent
\url{\codelink}

\section{Taskset Generation: Easy or Hard?}
\label{sMethod}

This section describes \ourmethod, a method for generating tasksets without requiring additional human knowledge, and which provides fine-grained control over the transfer difficulty between train and test tasks.

Our goal is to partition classes into train and test sets, such that we can control the difficulty of transferring a model trained on one set to the other.
At a high-level, \ourmethod finds those class partitions based on their \emph{class embedding} $\phi_1, \dots, \phi_M$.
Each class embedding $\phi_i$ is obtained by averaging the embedding of all samples from its corresponding class $y_i$:
\begin{equation*}
    \phi_i = \Exp[\mathbf{x} \sim p(\cdot \mid y = y_i)]{\phi(\mathbf{x})},
\end{equation*}
were $\phi$ is a pretrained feature extractor.
In practice, we find that pretraining $\phi$ on the dataset to be partitioned leads to satisfying results, but simpler approaches might also work. (\emph{e.g}\onedot} \def\Eg{\emph{E.g}\onedot an off-the-shelf model pretrained on ImageNet)

The key feature of \ourmethod is that train and test clusters are encouraged to stay at a \emph{prespecified distance of each other} in the class embedding space.
Intuitively, the farther and tighter the test set, the more difficult it will be for few-shot methods to discriminate between those classes.
On the other hand, a test task is easily solved if its classes are dissimilar (\emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot pooly clustered) and each class is fairly similar to classes seen during training (\emph{i.e}\onedot} \def\Ie{\emph{I.e}\onedot train and test clusters are adjacent).

To formalize this intuition, we pretend that the class embeddings were sampled from the mixture of two distributions $p_\text{train}(y)$ and $p_\text{test}(y)$.
We model $p_\text{train}$ and $p_\text{test}$ as multinomials such that their density for class $y_i$ depends on its distance to the distribution's centroid:
\begin{align*}
    p_\text{train}(y = y_i) &= \text{softmax} (\vert\vert \phi_i - \mu_\text{train} \vert \vert^2), \quad\text{and} \\
    p_\text{test}(y = y_i) &= \text{softmax} (\vert\vert \phi_i - \mu_\text{test} \vert \vert^2),
\end{align*}
where the centroids $\mu_\text{train}$ and $\mu_\text{test}$ are learnable parameters.
Maximizing the log-likelihood of the mixture $\frac1{2} (p_\text{train} + p_\text{test})$ is identical to a soft K-means objective if $p_\text{train}$ and $p_\text{test}$ were Gaussians~\cite[Chapter 22]{MacKay2002-la}.
In practice, we preferred multinomials because Gaussians led to numerical instabilities when optimizing the penalized objective (see below) by gradient descent.

We can incentivize $p_\text{train}$ and $p_\text{test}$ to lie at a given distance $R \geq 0$ by including a penalty $(D(p_\text{train} \vert\vert p_\text{test}) - R)^2$ based on our choice of statistical divergence $D$.
Combining this penalty with the mixture's log-likelihood yields our final penalized clustering objective:
\begin{align*}
    J = &- \sum_{i=0}^M \log \frac1{2} \left( p_\text{train}(y_i) + p_\text{test}(y_i)\right) \\
        &+ \lambda \left( D(p_\text{train} \vert \vert p_\text{test}) - R \right)^2,
\end{align*}
where $\lambda \geq 0$ balances the value of clustering versus the penalty.
Assuming train and test classes tightly cluster, controlling the transfer difficulty amounts to controlling the distance $D(p_\text{train} \vert\vert p_\text{test})$ between the two distributions, which we can achieve by adjusting $R$.

\paragraph{Train, Validation, and Test Assignments}
When $J$ is minimized with respect to $\mu_\text{train}$ and $\mu_\text{test}$, we get a solution we can use to partition $y_1, \dots, y_M$.
One approach is to assign class $y_i$ to the train set if the ratio $\frac{p_\text{train}(y_i)}{p_\text{test}(y_i)}$ is greater than 1, and to the test set if less.\footnote{Ties are broken at random, so that if $D(p_\text{train} \vert \vert p_\text{test}) = 0$ assignments are also random.}

Partitioning according to this decision rule is theoretically sound but might lead to slightly degenerate solutions in the context of few-shot learning.
For example, we observed instances where only 3 classes were assigned to the test set, which prevents evaluation in the 5-ways setting.
It also doesn’t prescribe how to devise a validation set from the test distribution.
To resolve those issues, we select the top 60\% scoring classes and assign them to the train set.
The remaining 40\% classes are split among validation and test sets in turn, according to their score: the least scoring class is assigned to the test set, the second-least scoring to the validation, etc...
Thus, validation and test classes retain roughly equal probability under $p_\text{test}$ and the resulting tasks are of similar difficulty.

\subsubsection*{Acknowledgements}
Fei Sha is on leave from University of Southern California.
This work is partially supported by NSF AwardsIIS-1513966/ 1632803/1833137, CCF-1139148, DARPA Award\#: FA8750-18-2-0117, DARPA-D3M - Award UCB-00009528, Google Research Awards, gifts from Facebook and Netflix, and ARO\# W911NF-12-1-0241 and W911NF-15-1-0484.

\clearpage
\pagebreak

{\small
\bibliographystyle{ieee_fullname}

\section{Related Works}
\label{sRelatedWorks}

\paragraph{Few-Shot Learning}
The goal of few-shot learning is to produce a model able to solve new tasks with access to only limited amounts of data~\cite{Fei-Fei2006-od, Miller2000-hy}.
It is closely related to meta-learning -- devising models that learn to learn~\cite{Bengio1991-yp, Schmidhuber1987-cx} -- but with a particular emphasis on the small quantity of available data.
This research direction has received a lot of interest in recent years, due to the numerous applications in natural language processing~\cite{Mu2019-hm,Yu2018-ns, Xu2020-jl}, medicine~\cite{Cai2020-as, Altae-Tran2017-uj, Prabhu2019-wi}, and more~\cite{Xing2019-ku, Bertinetto2016-cg, Brock2018-qy, Wei2019-hh}.
In the computer vision domain, a wide-range of approaches were proposed to tackle the few-shot image classification challenge.
Those include learning fast optimization schemes~\cite{Ravi2017-xu}, weight-imprinting~\cite{Qi2018-rw}, memory-augmented neural networks~\cite{Santoro2016-br}, casting the problem as stochastic process learning~\cite{Garnelo2018-pq}, or going about it from a Bayesian perspective\cite{Ortega2019-qi}.
This manuscript focuses on two major families of few-shot learning algorithms which we review below; for a more exhaustive account of few-shot and meta-learning, we refer the reader to the surveys of~\cite{Wang2020-fj} and~\cite{Vanschoren2018-mt}.

\paragraph{Gradient-Based Few-Shot Learning}
The main idea behind gradient-based few-shot learning algorithms is to discover a model whose parmeters can be adapted with just a few steps of gradient descent~\cite{Finn2018-fp}.
The representative algorithm of this family, MAML~\cite{Finn2017-gw}, does so by learning an initialization end-to-end while adaping all weights of the model.
Due to its flexible formulation, MAML was successfully applied to vision~\cite{Finn2017-gw}, robotics~\cite{Nagabandi2018-cz}, lifelong learning~\cite{Finn2019-yr}, and more~\cite{Zhou2020-ea}.
Variations of MAML improved upon its adaptation ability, reduce its computational footprint, or both.
Notable attempts to improve MAML's performance include meta-learning parameters dedicated to optimization, implicit~\cite{Lee2018-gw, Flennerhag2019-fw} or explicit~\cite{Li2017-li, Park2019-or, Behl2019-ki}, probabilistic regularization schemes~\cite{Grant2018-lh, Yoon2018-om}, as well as various training refinements~\cite{Antoniou2018-vl}.
To reduce the burden of computing second-order derivatives induced by the MAML objective, the authors suggested omitting those derivatives altogether at the cost of decreased performance.
When taking several gradient steps, one can leverage the implicit function theorem to more cheaply obtain the model updates~\cite{Rajeswaran2019-hx}.
Other options to mitigate the expense of second-order derivatives consists of conditioning the model on a latent embedding, and only updating that embedding during adaptation~\cite{Rusu2018-op, Rakelly2019-vk, Zintgraf2018-ot}.
In the same vein, the authors of ANIL~\cite{Raghu2019-ff} suggest that it is sufficient to update the very last layer of neural network to reap the benefits of MAML; that is, the model learns a feature extractor shared across all tasks and simply update the linear classifier a few times for each new task.

\paragraph{Metric-Based Few-Shot Learning}
Similar to ANIL, metric-based methods also share a feature extractor across tasks to extract high-dimensional embedding representations.
But, rather than adapting a linear classifier, they compute the distance of new instances' embeddings to the few-shot set of reference embeddings, akin to nearest neighbour classification.
Matching Networks~\cite{Vinyals2016-cm}, which was proposed to tackle the special case of one-shot learning, uses the (negative) cosine similarity as distance measure and relies on two different networks, one for the query and one for the reference embedding.
Prototypical Networks (ProtoNet)~\cite{Snell2017-kv} generalizes Matching Networks to the few-shot setting and measures similarity between embeddings with the Euclidean distance.
Some extensions to ProtoNet centered around learning and improving a specialized distance metric for a given task, typically by solving a convex optimization problem~\cite{Bertinetto2018-gn, Lee2019-cl, Zhang2020-ks}.
Others attempted to improve the embedding representations of a task~\cite{Oreshkin2018-yi, Rodriguez2020-vm, Ye2020-wf, Lichtenstein2020-py}.
Finally, ProtoMAML~\cite{Triantafillou2019-bj} combined MAML with ProtoNet such that the classification layer is initialized with the reference embedding values and the model is adapted for a few steps of gradient descent.

A recent line of work put a lot of those advances in question, and show that a simple baseline matches -- and often surpasses -- state-of-the-art performance in few-shot classification~\cite{Tian2020-gv, Wang2019-fc, Dhillon2019-mg, Chen2020-yd}.
Rather than learning by sampling many different tasks, they suggest to aggregate the classes from all tasks and to train a model via standard one-vs-rest multiclass classification.
At test time, the last linear layer of the model is removed to obtain a trained feature extractor, which can be used \`a la ProtoNet.
Directly~\cite{Tian2020-gv} or not, those results beg the question: \emph{do we still need to adapt our features for few-shot learning?}

\paragraph{Tasksets for Few-Shot Learning}
Our manuscript tackles this question when transfer to new tasks is especially challenging.
In this context, few-shot classification tasksets can be broadly categorized in two groups: those that leverage additional human information and those that don't.
In the latter case, the standard approach consists in taking an existing classification dataset and randomly partitioning classes in train and test sets.
Typical examples include CIFAR-FS~\cite{Bertinetto2018-gn}, mini-ImageNet~\cite{Vinyals2016-cm, Ravi2017-xu}, TCGA~\cite{Samiei2019-rx}, MultiMNIST~\cite{Sabour2017-ml, Sener2018-uc}, CU Birds~\cite{Welinder2010-cj}, FGVC Planes~\cite{Maji2013-ct} and Fungi~\cite{Sulc2020-vh, Schroeder2018-fu}, VGG Flowers~\cite{Nilsback2006-ye}, and Omniglot~\cite{Lake2015-uv, Lake2019-ph}.
Tasksets that do take advantage of human knowledge (usually in the form of semantic class relationships) attempt to minimize \emph{information overlap} to challenge transfer.
To the best of our knowledge, there are two such tasksets: FC100~\cite{Oreshkin2018-yi} and tiered-ImageNet~\cite{Ren2018-hk}.
The former is built on CIFAR100~\cite{Krizhevsky2009-mj} and leverages its superclass structure, while tiered-Imagenet is a subset of ImageNet-1k~\cite{Russakovsky2015-he} and uses the WordNet~\cite{Miller1990-zm} database.

Building on those datasets (and some more) is Meta-Dataset~\cite{Triantafillou2019-bj}: a large collection of tasksets aimed at mimicking transfer in real-world scenarios.
For example, one benchmarking scenario consists of training a model on ImageNet-1k and evaluating its performance on the remaining 9 tasksets.
Similar in spirit is VTAB~\cite{Zhai2019-wt}, which also consists of a large collection of synthetic and natural images aimed at evaluating representation learning and transfer.

In contrast to prior work, we study the question raised by Tian \emph{et al.}~\cite{Tian2020-gv}: is a good embedding sufficient for few-shot learning?
We approach this question from the standpoint where transfer is particularly challenging.
Since training gradient-based methods on large benchmarks such as Meta-Datasets and VTAB is prohibitively expensive, we devise a method to automatically generate tasksets with explicit control over transfer difficulty.


\section{Introduction}
The nonlinear Schr\"{o}dinger equation (NLSE) \cite{NLSE1,NLSE2,NLSE3,NLSE4,NLSE5,NLSE6,NLSE7,NLSE8,NLSE8b,NLSE8c}, which can be used to describe the dynamics of optical solitons in laser beam or  matter-wave solitons in Bose-Einstein condensate (BEC), has drawn more and more attention in recent decades. It is also widely used in quantum mechanics \cite{NLSE9}, hydrodynamics \cite{NLSE10}, nonlinear optics \cite{NLSE3,NLSE7}, BEC \cite{NLSE3,NLSE7} as well as superconductivity \cite{NLSE11}. Note that it is easy to stabilize a one-dimensional soliton, since it is now widely known that the bright soliton is the exact solution of one-dimensional NLSE with uniform cubic self-focusing nonlinearity. However, the collapse will happen in two-dimensional NLSE with uniform cubic self-focusing nonlinearity \cite{COLLAPSE}. A popular way to avoid such collapse is the introduction of the linear potentials \cite{LP1,LP2,LP3}, which can stabilize various kinds of solitons in all dimensions. Another way to create stable solitons in high dimensions is the introduction of nonlinearities, such as the employment of periodically modulated nonlinearities \cite{NLSE3,MN1,MN2} or spatially inhomogeneous defocusing nonlinearities \cite{DEF1,DEF2,DEF3,DEF4,DEF5}.

It is worth to mention that single NLSE can only describe the propagation of an optical soliton in a laser beam or the evolution of a matter-wave soliton in a BEC. As for the settings of two (or more) laser beams or BECs, the models of coupled NLSEs \cite{CNLSE1,CNLSE2,CNLSE3,CNLSE4,CNLSE5,CNLSE6,CNLSE7} have to be considered due to the coupled interactions of the laser beams or BECs. These models are widely studied in optical couplers, especially in the fields of optical fiber \cite{CNLSE8,CNLSE9} and waveguide \cite{CNLSE10} in recent years. Further, the models of couplers with parity-time- ($\mathcal{PT}$-) symmetry have drawn more and more attention, since Bender demonstrated that non-Hermitian Hamiltonians with $\mathcal{PT}$-symmetry can have real spectra \cite{PT0}, which is an important extension of standard quantum mechanics. These physical models with $\mathcal{PT}$-symmetry exhibit some novel physical properties due to their non-Hermitian Hamiltonians \cite{PT1,PT2} and various physical settings with $\mathcal{PT}$-symmetry have been widely investigated during the past decade \cite{PT3,PT4,PT5,PT6,PT7,PT8,PT9}.

The fractional derivative \cite{Frac}, firstly introduced by Laskin, is another great extension of the standard quantum mechanics \cite{Lask1,Lask2,Lask3}. Under such condition, the quantum mechanics becomes the fractional quantum mechanics, and the NLSE becomes the nonlinear fractional Schr\"{o}dinger equation (NLFSE). Followed by the above landmark study of Laskin in the NLSE, various studies have been reported in recent years \cite{Frac00,Frac0,Frac1,Frac2,Frac3,Frac4,Frac5,Frac6,Frac7,Frac8,Frac9,Frac10,Frac11,Frac12,Frac13,Frac14,Frac15,Frac15b,Frac16,Frac17,Frac18,Frac19,Frac19b,Frac19c,Frac20,Frac21,Frac22,Frac23,Frac24,Frac25,Frac26,Frac27,Frac28}, including the propagation of light beams \cite{Frac1,Frac2}, accessible solitons \cite{Frac5,Frac6}, gap solitons \cite{Frac7,Frac16}, vortex solitons \cite{Frac18,Frac19c,Frac27}, solitons in nonlinear lattices \cite{Frac14}, and soliton clusters \cite{Frac19,Frac19b}.

Despite so many researches that have been reported in the NLFSE, the coupled $\mathcal{PT}$-symmetric NLFSEs with gain and loss, which can be used to describe the propagations of binary $\mathcal{PT}$-symmetric laser beams or BECs with coupled interactions, have not been yet reported, to the best of our knowledge. This coupled model can be employed to study the optical waveguides and fiber couplers, which have great potential in optical communications. The objective of this work is to extend the NLFSE to its coupled form with $\mathcal{PT}$-symmetry, from which we can explore the existence, stability and other properties of optical and matter-wave solitons in this model under fractional-order diffraction.

This paper is organized as follows. In Sec. \ref{sec2}, we introduce the theoretical model and report some typical analytical solutions for coupled solitons. Sec. \ref{sec3}, which presents numerical results for $\mathcal{PT}$-symmetric coupled solitons, is divided into two parts. The numerical results for the model without cross-interactions modulations are given in Sec. \ref{sec3a}, and the corresponding numerical results for the model with cross-interactions modulations are reported in Sec. \ref{sec3b}. We find that the field profiles and stability domains of coupled solitons are strongly affected by the linear coupling coefficient, the fractional-order diffraction described by the L\'{e}vy index $\alpha$, the propagation constant, and the cross-interactions modulation. It should be mentioned that these coupled solitons can exist only when the gain/loss parameter ($\gamma$) is less than the linear coupling coefficient ($\kappa$). Finally, the paper is summarized in Sec. \ref{sec4}.

\section{Theoretical model}
\label{sec2}
The propagations of laser beams in $\mathcal{PT}$-symmetric couplers with gain and loss under fractional-order diffraction can be described by the dimensionless coupled NLFSEs
\begin{equation}
\left\{
\begin{aligned}
i\frac{\partial U_1}{\partial z}=&\frac{1}{2}\left(-\frac{\partial^2}{\partial x^2}\right)^{\alpha/2}U_1-\mathrm{g}\left|U_1\right|^2U_1-\epsilon\left|U_2\right|^2U_1 \\
&-\kappa U_2+i\gamma U_1, \\
i\frac{\partial U_2}{\partial z}=&\frac{1}{2}\left(-\frac{\partial^2}{\partial x^2}\right)^{\alpha/2}U_2-\mathrm{g}\left|U_2\right|^2U_2-\epsilon\left|U_1\right|^2U_2 \\
&-\kappa U_1-i\gamma U_2.
\end{aligned}
\label{NLFSE}
\right.
\end{equation}
Here $U_{1,2}$ and $z$ represent the field amplitudes and propagation distance, respectively, $\mathrm{g}>0$ and $\epsilon>0$ represent the coefficients of attractive self- and cross-interactions modulations of the components, respectively. Note that $\mathrm{g}<0$ and $\epsilon<0$ denote the coefficients of repulsive self- and cross-interactions modulations of the components that can be used to generate dark solitons \cite{DARK}. In this article, we  focus on the formation of bright solitons, hence the case when $\mathrm{g}<0$ and $\epsilon<0$ is not discussed here. The parameter $\kappa>0$ denotes the linear coupling  coefficient, and $\gamma>0$, that is the coefficient of gain and loss, stands for the gain in one component and the loss in the other one.   Note that $(-\partial^2/\partial x^2)^{\alpha/2}$ denotes the fractional derivative, where $\alpha$ ($1<\alpha\leq2$) stands for the L\'{e}vy index. The definition of fractional derivative is as follows \cite{Lask1,Lask2,Lask3}
\begin{equation}
\begin{split}
\left( -\frac{\partial ^{2}}{\partial x^{2}}\right) ^{\alpha /2}U=&\frac{1}{2\pi }\int_{-\infty }^{+\infty }|s|^{\alpha }ds \\
\times & \int_{-\infty }^{+\infty}d\zeta \exp \left[ ip\left( x-\zeta \right) \right] U(\zeta ).
\end{split}
\end{equation}
In particular, Eqs. (\ref{NLFSE}) will degenerate to the generic $\mathcal{PT}$-symmetric coupled nonlinear Schr\"{o}dinger equations when $\alpha=2$. As for BECs, $U_{1,2}$ and $z$ should be replaced by the wave functions $\phi_{1,2}$ and time $t$.

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig1.eps}
\end{center}
\caption{Profiles of solitons (shown for components of $U_1$): (a) with different values of $\alpha$ at $\kappa=0.8$, $\gamma=0.001$, $b=1.2$; (b) with different values of $b$ at $\alpha=1.8$, $\kappa=0.8$, $\gamma=0.001$; (c) with different values of $\kappa$ at $\alpha=1.8$, $\gamma=0.001$, $b=1$; (d) with different values of $\gamma$ at $\alpha=1.8$, $\kappa=1$, $b=1.5$. The profiles of $|U_2|$ are similar to $|U_1|$. $\epsilon=0$ are used in Figs. \ref{fig1}$\sim$\ref{fig3}. $\mathrm{g}=1$ are used throughout this paper.}
\label{fig1}
\end{figure}

With the real propagation constant $b$ ($b$ should be replaced by chemical potential $-\mu$ in BECs), the stationary solutions of these $\mathcal{PT}$-symmetric coupled solitons with gain and loss can be found by generic forms
\begin{equation}
\left\{
\begin{aligned}
U_1(x,z)=U(x){\rm exp}(ibz-i\delta/2),\\
U_2(x,z)=U(x){\rm exp}(ibz+i\delta/2).
\end{aligned}
\label{NLFSES}
\right.
\end{equation}
Here $\delta$ denotes the constant phase shift between the components, which is given by
\begin{equation}
\delta={\rm arcsin}(\gamma/\kappa).
\label{CPS}
\end{equation}
As is seen from Eqs. (\ref{NLFSES}) and (\ref{CPS}), the difference between the two components $U_{1,2}$, which is caused by their parameter of gain-loss, is their imaginary part ${\rm exp}(\pm i\delta/2)$. Due to their opposite parameter of gain-loss, these two components will have the opposite transition of energy in their propagations, which can be referred to the following Figs. \ref{fig3}, \ref{fig6}, \ref{fig7}, \ref{fig8} and \ref{fig9}.

Here the real function $U(x)$  satisfies the generic stationary NLFSE with a varying value of the propagation constant $b$
\begin{equation}
-(b-\kappa_s)U=\frac{1}{2}\left(-\frac{\partial^2}{\partial x^2}\right)^{\alpha/2}U-\mathrm{g} U^3-\epsilon U^3,
\label{NLFSE2}
\end{equation}
\begin{equation}
\kappa_s=\sqrt{\kappa^2-\gamma^2}.
\label{KS}
\end{equation}
According to Eq. (\ref{KS}), it is easy to see that these coupled solitons can exist only when
\begin{equation}
\gamma\leq\gamma_{max}=\kappa.
\label{GAMMAK}
\end{equation}
From Eq. (\ref{GAMMAK}), the largest value of the gain-loss coefficient in these $\mathcal{PT}$-symmetric couplers is limited to the value of the linear coupling coefficient $\kappa$.

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig2.eps}
\end{center}
\caption{(a) Soliton power $P$ versus $\alpha$ at $\kappa=0.8$, $\gamma=0.001$, $b=1.2$. (b) $P$ versus $b$ at $\alpha=1.8$, $\kappa=0.8$, $\gamma=0.001$. (c) $P$ versus $\kappa$ at $\alpha=1.8$, $\gamma=0.001$, $b=1$. (d) Stability (white) and instability (gray) domains for the coupled solitons with $\kappa=0.8$, $\gamma=0.001$ in the ($b,\alpha$) plane. The propagations of the solitons corresponding to the points marked by A1 and A2 are shown in Figs. \ref{fig3}(a,b) and \ref{fig3}(c,d), respectively.}
\label{fig2}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig3.eps}
\end{center}
\caption{Propagations of coupled solitons with different values of $b$ at $\alpha=1.8$, $\kappa=0.8$, $\gamma=0.001$: (a,b) with $b=1$; (c,d) with $b=1.35$. The left~/~right column displays the $U_1$~/~$U_2$ component.}
\label{fig3}
\end{figure}

Particularly, analytical solution of Eq. (\ref{NLFSE2}) can be solved when $\alpha=2$, which is given by
\begin{equation}
U=\sqrt{2(b-\kappa_s)/(\mathrm{g}+\epsilon)}{\rm sech}\left(\sqrt{2(b-\kappa_s)}x\right),
\label{SOLUTION}
\end{equation}
where $\kappa_s$ is defined by Eq. (\ref{KS}).

The soliton power $P$ of $U_{1,2}$ is defined by $P_1=\int\left|U_1\right|^2dx$, $P_2=\int\left|U_2\right|^2dx$. In this paper, the modified squared-operator method \cite{MSOM} and split-step Fourier method are employed to solve the stationary solutions for Eq. (\ref{NLFSE2}) and numerical simulations for Eqs. (\ref{NLFSE}) respectively.

\section{Numerical results}
\label{sec3}

In this section, we focus on presenting our numerical results of $\mathcal{PT}$-symmetric coupled solitons with gain and loss under different parameters in the models of Eqs. (\ref{NLFSE}) and Eq. (\ref{NLFSE2}), inlcuding the $\mathcal{PT}$-symmetric coupled solitons with different linear coupling coefficients, fractional-order diffractions, propagation constants and cross-interactions modulations. Particularly, we find that the cross-interactions modulation can greatly affect these coupled solitons. Thus we divide our numerical results into two parts, namely the coupled solitons without/with cross-interactions modulations which can be referred to Sec. \ref{sec3a} and Sec. \ref{sec3b}, respectively.

\subsection{Coupled solitons without cross-interactions modulation}
\label{sec3a}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig4.eps}
\end{center}
\caption{Profiles of solitons (shown for components of $U_1$): (a) with different values of $\alpha$ at $\kappa=0.8$, $\epsilon=0.5$, $\gamma=0.001$, $b=2$; (b) with different values of $b$ at $\alpha=1.7$, $\kappa=0.8$, $\epsilon=0.5$, $\gamma=0.001$; (c) with different values of $\kappa$ at $\alpha=1.7$, $\epsilon=0.5$, $\gamma=0.001$, $b=2$; (d) with different values of $\epsilon$ at $\alpha=1.7$, $\kappa=0.8$, $\gamma=0.001$, $b=2$. The profiles of $|U_2|$ are similar to $|U_1|$.}
\label{fig4}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig5.eps}
\end{center}
\caption{(a) $P$ versus $\alpha$ at $\kappa=0.8$, $\epsilon=0.5$, $\gamma=0.001$, $b=2$. (b) $P$ versus $b$ at $\alpha=1.7$, $\kappa=0.8$, $\epsilon=0.5$, $\gamma=0.001$. (c) Stability (white) and instability (gray) domains for the coupled solitons with $\alpha=1.7$, $\gamma=0.001$, $b=1.5$ in the ($\epsilon,\kappa$) plane. The propagations of the solitons corresponding to the points marked by B1 and B2 are shown in Figs. \ref{fig6}(a,b) and \ref{fig6}(c,d), respectively. (d) Stability (white) and instability (gray) domains for the coupled solitons with $\kappa=0.8$, $\epsilon=0.5$, $\gamma=0.001$ in the ($b,\alpha$) plane. The propagations of the solitons corresponding to the points marked by C1 and C2 are displayed in Figs. \ref{fig7}(a,b) and \ref{fig7}(c,d), respectively.}
\label{fig5}
\end{figure}

The $\mathcal{PT}$-symmetric physical model without cross-interactions modulation means that the coefficient $\epsilon=0$, which is a simple model but a nontrivial one. We now display the relevant numerical results in this setting. Note that the numerical results under this setting are shown in Figs. \ref{fig1}$\sim$\ref{fig3}. According to Eq. (\ref{KS}), it is not difficult to find that these real $\mathcal{PT}$-symmetric coupled solutions can only be supported under the condition of $\gamma\leq\kappa$. For simplicity, note also that we use $\mathrm{g}=1$ throughout this article. Profiles of the $\mathcal{PT}$-symmetric coupled solitons with different values of the L\'{e}vy index $\alpha$, different propagation constants  $b$, different linear coupling coefficients $\kappa$, and different values of the gain and loss parameter $\gamma$ are shown in Figs. \ref{fig1}(a$\sim$d), respectively. According to Eqs. (\ref{NLFSES}), we see that the stationary solutions of solitons $U_1(x)$ and $U_2(x)$ have the same real part $U(x)$ and conjugate imaginary parts ${\rm exp}(\pm i\delta/2)$, that is, $U_1(x)=U(x){\rm exp}(-i\delta/2)$, $U_2(x)=U(x){\rm exp}(i\delta/2)$. Note also that the propagation constant $b$ will disappear in the expressions of the stationary solutions of solitons $U_1(x)$ and $U_2(x)$ after we submit Eqs. (\ref{NLFSES}) into Eqs. (\ref{NLFSE}). According to the above expressions of $U_1(x)$ and $U_2(x)$, the profiles of the modulus $|U_2|$ are similar to those of the modulus $|U_1|$, thus we display here only the profiles of $|U_1|$. From Fig. \ref{fig1}(a), the amplitude of coupled solitons decreases and the width increases when $\alpha$ increases. According to Fig. \ref{fig1}(b), the amplitude of coupled solitons increases and the width decreases when $b$ increases. Note also that the amplitude of coupled solitons decreases and the width increases when $\kappa$ increases, as presented in Fig. \ref{fig1}(c), similar to the situation in Fig. \ref{fig1}(a). In Fig. \ref{fig1}(d), the amplitude of coupled solitons increases with the increase of $\gamma$ and the width of solitons decreases with the increase of $\gamma$.

We now focus on the soliton power of the obtained families of coupled solitons, including the relationships of the soliton power $P$ versus the L\'{e}vy index $\alpha$, $P$ versus the propagation constant $b$, and $P$ versus the linear coupling coefficient $\kappa$. We point out that the dependences of the soliton power of the $|U_2|$ component on the  parameters of the model are similar to those of the $|U_1|$ component. Fig. \ref{fig2}(a) shows that the soliton power $P_1$ increases with the increase of $\alpha$, and the increasing speed gradually decreases. According to Fig. \ref{fig2}(b), the soliton power increases when $b$ increases, which shows that the relationship between soliton power and propagation constant satisfies the Vakhitov-Kolokolov (VK) stability criterion $\partial P/\partial b >0$, which is a necessary but not sufficient condition for stable soliton solutions \cite{VK}. The soliton power decreases with the increase of $\kappa$, as presented in Fig. \ref{fig2}(c). The stable and unstable domains of these $\mathcal{PT}$-symmetric coupled solitons are very important, here, Fig. \ref{fig2}(d) displays the stability (white) and instability (gray) domains for the obtained coupled solitons in the ($b,\alpha$) plane. From Fig. \ref{fig2}(d), we can find that the coupled solitons can be stable only when the L\'{e}vy index $\alpha\geq1.2$. The propagation dynamics of the coupled solitons marked by the points A1 and A2 are presented in Figs. \ref{fig3}(a,b) and \ref{fig3}(c,d), respectively. The coupled solitons marked by the point A1 are stable, and both components of the coupled solitons are keeeping their amplitudes and shapes even after a long propagation distance (propagation distance $z=1000$), as displayed in Figs. \ref{fig3}(a,b). On the other hand, the coupled solitons marked by the point A2 are unstable, and both components of the coupled solitons are greately distorted after some propagation distance, as presented in Figs. \ref{fig3}(c,d).

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig6.eps}
\end{center}
\caption{Propagations of coupled solitons with different values of $\kappa$ at $\alpha=1.7$, $\epsilon=0.7$, $\gamma=0.001$, $b=1.5$: (a,b) with $\kappa=0.8$; (c,d) with $\kappa=0.2$. The left~/~right column displays the $U_1$~/~$U_2$ component.}
\label{fig6}
\end{figure}

\subsection{Coupled solitons with cross-interactions modulation}
\label{sec3b}

In this subsection, we focus on coupled solitons with cross-interactions modulation, that is, when $\epsilon>0$). It should be mentioned that the numerical results under this setting are displayed in Figs. \ref{fig4}$\sim$\ref{fig9}. Fig. \ref{fig4} shows the profiles of $\mathcal{PT}$-symmetric coupled solitons with different parameters, including the solitons with different L\'{e}vy indices $\alpha$, different propagation constants $b$, different linear coupling coefficients $\kappa$, and different cross-interactions modulations $\epsilon$. The profiles of $|U_2|$ are similar to $|U_1|$, therefore, we only display the profiles of $|U_1|$ in this subsection.

The amplitude of coupled solitons decreases and the width of increases when $\alpha$ increases, as shown in Fig. \ref{fig4}(a). In Fig. \ref{fig4}(b), the amplitude of coupled solitons increases while the width decreases with the increase of $b$. According to Fig. \ref{fig4}(c), the amplitude of coupled solitons decreases and the width increases when $\kappa$ increases. From Fig. \ref{fig4}(d), the amplitude and width of coupled solitons increase when $\epsilon$ increases. Obviously, the results shown in Figs. \ref{fig4}(a$\sim$c) are similar to their counterparts shown in Figs. \ref{fig1}(a$\sim$c) (without cross-interactions modulation).

Fig. \ref{fig5} displays the soliton power and relevant stability area of these $\mathcal{PT}$-symmetric coupled solitons with cross-interactions modulation. Here, the soliton power of $|U_2|$ is similar to that of $|U_1|$. Fig. \ref{fig5}(a) demonstrates that the soliton power $P_1$ increases with the increase of $\alpha$, and the growth speed gradually diminishes. From Fig. \ref{fig5}(b), the soliton power increases when $b$ increases, which demonstrates that the relationship between soliton power and propagation constant for the obtained coupled solitons also satisfies the VK stability criterion. The stable and unstable domains of the coupled solitons with cross-interactions modulation are also very important. Here, Figs. \ref{fig5}(c,d) show the stability (white) and instability (gray) domains for the coupled solitons in the ($\epsilon,\kappa$) and ($b,\alpha$) planes, respectively. According to Fig. \ref{fig5}(c), the coupled solitons can be stable only when $\kappa$ is above a certain threshold, which is about 0.55 in this case. Fig. \ref{fig5}(d) demonstrates that these coupled solitons can be stable only when $\alpha\geq1.2$. Further, by comparing the results of Fig. \ref{fig5}(d) and Fig. \ref{fig2}(d), we can clearly see that the stable domain of Fig. \ref{fig5}(d) is much broader than the one in Fig. \ref{fig2}(d). Note that other parameters in Fig. \ref{fig5}(d) and Fig. \ref{fig2}(d) are the same except for cross-interactions modulation $\epsilon$, therefore, we can draw a conclusion that $\epsilon$ can greatly affect the stability and instability domains.

The propagation dynamics of the coupled solitons marked by the points (B1,B2) [in Fig. \ref{fig5}(c)] and (C1,C2) [in Fig. \ref{fig5}(d)] are presented in Figs. \ref{fig6} and \ref{fig7}, respectively. The coupled solitons marked by the points B1 and C1 are stable, and both components of these coupled solitons are keeping their amplitudes and shapes even after long-distance propagations (after a propagation distance $z=1000$), as displayed in Figs. \ref{fig6}(a,b) and \ref{fig7}(a,b). On the other hand, the coupled solitons marked by the points B2 and C2 are unstable, and both components of these coupled solitons are greately distorted even after short-distance propagations, as presented in Figs. \ref{fig6}(c,d) and \ref{fig7}(c,d).

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig7.eps}
\end{center}
\caption{Propagations of coupled solitons with different values of $b$ at $\alpha=1.6$, $\kappa=0.8$, $\epsilon=0.5$, $\gamma=0.001$: (a,b) with $b=1.5$; (c,d) with $b=2.4$. The left~/~right column displays the $U_1$~/~$U_2$ component.}
\label{fig7}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig8.eps}
\end{center}
\caption{Stable propagations of oscillating solitons at $\alpha=1.8$, $\epsilon=1$, $b=2$: (a,b) with $\kappa=0.02$, $\gamma=0.01$; (c,d) with $\kappa=0.04$, $\gamma=0.01$; (e,f) with $\kappa=0.04$, $\gamma=0.038$. The left~/~right column displays the $U_1$~/~$U_2$ component. Note that the oscillating period can be controlled only by varying the value of $\kappa$. The soliton power $P$ versus propagation distant $z$ for panels (a)$\sim$(c) are displayed in Figs. \ref{fig9}(a)$\sim$\ref{fig9}(c), respectively.}
\label{fig8}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\includegraphics[width=1\columnwidth]{Fig9.eps}
\end{center}
\caption{Soliton power $P$ versus $z$ for the propagation displayed in Figs. \ref{fig8}(a)$\sim$\ref{fig8}(c). The yellow and blue lines stand for $P_1$ and $P_2$, respectively. Note that the parameter $\gamma$ can not affect the oscillating period ${\rm T}$, while it can affect the convergence speed to the equilibrium state.}
\label{fig9}
\end{figure}

Next we turn to investigate the formation and propagation properties of the oscillating $\mathcal{PT}$-symmetric coupled solitons in the model of cross-interactions modulation. It was mentioned before that $\mathrm{g}=1$ is used throughout this paper; here we find that the unstable coupled solitons become oscillating coupled solitons when $\epsilon=1$. Some examples of the oscillating coupled solitons are displayed in Fig. \ref{fig8}, which shows the propagation dynamics of $U_1$ in the left column and the corresponding propagation dynamics of $U_2$ in the right column. Figs. \ref{fig8}(a,b) display the propagation dynamics of the oscillating coupled solitons with a long oscillating period, and Figs. \ref{fig8}(c$\sim$f) present the propagation dynamics of the oscillating coupled solitons with a short oscillating period. According to Figs. \ref{fig8}(a,b), we see that when the amplitude of one component decreases, the amplitude of the other component would increase at the same time, which is related to the coefficient of gain and loss $\gamma$.

To further study the properties of these interesting oscillating $\mathcal{PT}$-symmetric coupled solitons, we draw the relationships between their soliton power $P$ and propagation distance $z$ in Fig. \ref{fig9} (whose propagation dynamics are presented in Fig. \ref{fig8}). In Fig. \ref{fig9}, the yellow and blue lines stand for P1 and P2, respectively. By comparing the results in Figs. \ref{fig9}(a,b), we can see that the oscillating period depends on the linear coupling coefficient $\kappa$, and it decreases when $\kappa$ increases. In addition, we point out that the convergence speed of transforming the coupled solitons into oscillating coupled solitons is strongly related to the specific value of the gain/loss parameter $\gamma$. That is to say, the larger $\gamma$ is, the faster the coupled solitons transform into the oscillating coupled solitons. It should be emphasized that we have verified the above conclusions by using different sets of parameters. Furthermore, we find that these oscillating $\mathcal{PT}$-symmetric coupled solitons can be supported under the condition of $\mathrm{g}=\epsilon>0$, which is not limited to $\mathrm{g}=\epsilon=1$.

\section{Conclusion}
\label{sec4}
In this paper, we have investigated the existence and propagation dynamics of the families of coupled solitons in $\mathcal{PT}$-symmetric physical models with gain and loss in fractional dimension, including the coupled models with and without cross-interactions modulation. The profiles, powers, stability domains, and propagation properties of these $\mathcal{PT}$-symmetric coupled solitons are studied in both physical settings. According to our study, the profiles of such $\mathcal{PT}$-symmetric coupled solitons are sensitive to the values of propagation constant, and the stable coupled solitons can only exist under the condition that the L\'{e}vy index is greater than a critical value. In addition, the obtained results  reveal that the cross-interactions modulation can greatly affect the stability and instability domains of the coupled solitons. In other words, the stability area of the model with cross-interactions modulation is much broader than the one without cross-interactions modulation. Further, we have also studied the oscillating $\mathcal{PT}$-symmetric coupled solitons that can be supported when the coefficients of the self- and cross-interactions modulations are the same. Besides, the oscillating period of these coupled solitons can be controlled by varying the value of the linear coupling coefficient, and the convergence speed of coupled solitons that are transforming into oscillating coupled solitons is only related to the value of the gain/loss parameter.

The physical model studied in this work can be extended to the (2+1)-dimensional setting, which is an interesting but a quite difficult problem due to the presence of the critical collapse \cite{COLLAPSE}. Also, the introduction of linear periodic potentials into the physical model is another interesting issue to be investigated, because the coupled gap solitons \cite{CGS} can form in that setting.

\section*{Acknowledgment}
We thank Dumitru Mihalache and Jiawei Li for their help on improving the writing of the manuscript.

\section*{Funding}
This work was supported by National Major Instruments and Equipment Development Project of National Natural Science Foundation of China (No. 61827815), by National Natural Science Foundation of China (No. 62075138), and by Science and Technology Project of Shenzhen (Nos. JCYJ20190808121817100, JCYJ20190808164007485).

\section*{Conflict of interest}
\label{sec5}
The authors declare that they have no conflict of interest.


\section{Introduction}

For any $n \geq 1$, a partition of $[n]=\{1,\ldots,n\}$ is a collection of disjoint nonempty subsets of $[n]$ whose union is $[n]$.
The most common method of representing a set partition given by $\{B_1,B_2,\ldots,B_k\}$ is to write the blocks as
\begin{equation*}
    B_1/B_2/\cdots/B_k.
\end{equation*}
The usual convention followed is to order the blocks such that
\begin{equation*}
    \operatorname{min}(B_1) < \operatorname{min}(B_2) < \cdots < \operatorname{min}(B_k).
\end{equation*}Sometimes it is also convenient, when it does not cause confusion, to write the elements of each block without braces or commas.
For example, the partition $\{\{1,3,4\},\{2,6\},\{5,7,8\}\}$ of $[8]$ is written as $\{1,3,4\}/\{2,6\}/\{5,7,8\}$ or $134/26/578$.

A partition $B_1/B_2/\cdots/B_k$ of $[n]$ such that $\operatorname{min}(B_1) < \operatorname{min}(B_2) < \cdots < \operatorname{min}(B_k)$ can be represented by the sequence $a_1a_2\cdots a_n$ where $a_i=j$ if $i \in B_j$.
This sequence is called the \textit{restricted growth function} associated to the partition. For example, the partition $134/26/578$ has corresponding sequence $12113233$. For more on restricted growth functions, the interested reader is referred to Sagan's paper \cite[Section 4]{sag}. Henceforth, we will represent partitions using restricted growth functions.
The following definition of pattern avoidance in set partitions was first introduced by Sagan \cite{sag}.

\begin{definition}
A partition $\sigma=\sigma_1\sigma_2\cdots \sigma_n$ \textit{contains} a partition (or pattern) $\pi=\pi_1\pi_2\cdots\pi_m$ if there exists a subsequence $1 \leq h(1) < h(2) < \cdots < h(m) \leq n$ such that for any $i,j \in [m]$, $\sigma_{h(i)}=\sigma_{h(j)}$ if and only if $\pi_i=\pi_j$ and $\sigma_{h(i)}<\sigma_{h(j)}$ if and only if $\pi_i<\pi_j$. If $\sigma$ does not contain $\pi$, we say that $\sigma$ \textit{avoids} the pattern $\pi$.
\end{definition}

The topic of pattern avoidance has been an active area of research in enumerative combinatorics, starting with Knuth's work on permutations \cite{knu}. The study of pattern avoidance in set partitions was initiated by Klazar \cite{kla2}. Since then, several different notions of pattern avoidance of set partitions have been studied (see, e.g., the work of Chen et al. \cite{cddsy}, Goyt \cite{goyt}, or Bloom and Elizalde \cite{bloeli}).

Let $T$ be a set of patterns. We will denote the set of partitions of $[n]$ that avoid all the patterns of $T$ as $P_n(T)$ and the number of such partitions, i.e., $|P_n(T)|$ as $p_n(T)$.

\begin{definition}
Two sets of patterns $T$ and $R$ are said to be Wilf-equivalent, written as $T \sim R$, if $p_n(T)=p_n(R)$ for all $n \geq 1$.
\end{definition}

For example, $P_n(\{123\})$ is the set of partitions of $[n]$ with at most two blocks and $P_n(\{122\})$ is the set of partitions where any block not containing $1$ is a singleton.
Hence we get $\{123\} \sim \{122\}$ since $p_n(\{123\})=p_n(\{122\})=2^{n-1}$ for all $n \geq 1$.

In this article, our focus will be on pairs of patterns $\{\sigma,\tau\}$, where $\sigma$ is a pattern of size $3$ with at least $2$ blocks.
Such pairs are called $(3,k)$-pairs when the the size of $\tau$ is $k$.
Jel\'{\i}nek, Mansour, and Shattuck \cite{jel} studied such pairs and obtained an upper bound for the number of Wilf-equivalence classes of $(3,k)$-pairs.
This was done by describing various Wilf-equivalences among $(3,k)$-pairs.
They also showed that any other Wilf-equivalences, if they exist, are between $(3,k)$-pairs of the form $\{112,\tau\}$ and left open the following problem.

\begin{problem*}[{\cite[Problem 2.17]{jel}}]
Are there any more equivalences among the $(3,k)$-pairs of the
form $\{112,\tau\}$ other than those that we know about? Equivalently, are there any two distinct 2-free integer partitions that are $\overset{d}{\sim}$-equivalent (see \cref{defd})?
\end{problem*}

The main aim of this article is to answer this question.
In \cref{sec2}, we describe the results obtained by Jel\'{\i}nek, Mansour, and Shattuck \cite{jel} , and show how the second question is equivalent to the first in the above problem.
In \cref{sec3}, we answer this equivalent question and hence prove that there are no other Wilf-equivalences among $(3,k)$-pairs other than the ones described by Jel\'{\i}nek, Mansour, and Shattuck in \cite{jel}.

\section{Wilf-equivalences among $(3,k)$-pairs}\label{sec2}

In this section we will state the Wilf-equivalences between $(3,k)$-pairs derived by Jel\'{\i}nek, Mansour, and Shattuck \cite{jel}.

Note that if we are studying $P_n(T)$ where $T=\{\tau_1,\tau_2,\ldots\}$, we can assume that $\tau_i$ avoids $\tau_j$ for all $i \neq j$.
For example, when avoiding a $(3,k)$-pair $\{121,\tau\}$, we can assume that $\tau$ avoids $121$.
It can be shown that if the number of blocks and size of each block is specified, there is a unique partition that avoids $121$.
Namely, if a partition has $m$ blocks such that the $i^{th}$ block has $a_i$ elements for all $i \in [m]$, the unique such partition avoiding $121$ is $1^{a_1}2^{a_2}\cdots m^{a_m}$.
Here, $k^a$ represents $a$ consecutive copies of $k$.

\begin{definition}
A \textit{composition} of a positive integer $n$ is a sequence $(a_1,a_2,\ldots,a_m)$ such that $a_1+\cdots+a_m=n$.
The set of compositions of $n$ is denoted as $C_n$.
\end{definition}

For example, $C_3=\{(1,1,1),(1,2),(2,1),(3)\}$.
The above discussion implies that $P_n(\{121\})$ is in bijection with $C_n$.
The partition avoiding $121$ associated to the composition $a=(a_1,\ldots,a_m)$ is
\begin{equation*}
    \tau_{121}(a)=1^{a_1}2^{a_2}\cdots m^{a_m}.
\end{equation*}
It can be shown that the partitions of the form $\tau_{121}(b)$ that contain $\tau_{121}(a)$ are those where $b$ has a subsequence, having same length as $a$, such that each term in the subsequence has value at least that of the corresponding term in $a$.

\begin{definition}
The composition $b=(b_1,\ldots,b_k)$ is said to \textit{dominate} the composition $a=(a_1,\ldots,a_m)$ if there exists a subsequence $1\leq i(1) < i(2) < \cdots < i(m) \leq k$ such that $b_{i(j)} \geq a_j$ for all $j \in [m]$.
Such a subsequence is called an \textit{occurrence} of $a$ in $b$.
For any positive integer $n$, $D_n(a)$ will denote the set of compositions of $n$ that dominate $a$.
\end{definition}

Hence, we get that $P_n(\{121,\tau_{121}(a)\})=\{\tau_{121}(b) \mid b \in C_n \setminus D_n(a)\}$.
We now descibe the Wilf-equivalences among $(3,k)$-pairs of the form $\{121,\tau_{121}(a)\}$.
From the description of $P_n(\{121,\tau_{121}(a)\})$, it is clear that $\{121,\tau_{121}(a)\} \sim \{121,\tau_{121}(a')\}$ if and only if $|D_n(a)|=|D_n(a')|$ for all $n \geq 1$.
Hence we make the following definition.

\begin{definition}\label{defd}
For two compositions $a$ and $a'$, we say $a$ and $a'$ are \textit{dominating equivalent}, written as $a \overset{d}{\sim} a'$, if for all positive integers $n$, $|D_n(a)|=|D_n(a')|$.
\end{definition}

Therefore, $\{121,\tau_{121}(a)\} \sim \{121,\tau_{121}(a')\}$ if and only if $a \overset{d}{\sim} a'$.
Using similar ideas it can also be shown that $P_n(\{112\})$ is in bijection with compositions of $n$, where to a composition $a=(a_1,\ldots,a_m)$ we associate the partition
\begin{equation*}
    \tau_{112}(a)=12\cdots(m-1)mm^{a_m-1}\cdots2^{a_2-1}1^{a_1-1}.
\end{equation*}
Also, $P_n(\{112,\tau_{112}(a)\}) = \{\tau_{112}(b) \mid b \in C_n \setminus D_n(a)\}$ and $\{112,\tau_{112}(a)\} \sim \{112,\tau_{112}(a')\}$ if and only if $a \overset{d}{\sim} a'$.

Using that fact that $p_k(\sigma)=2^{k-1}$ for all patterns $\sigma$ of size $3$ other than $111$ (see \cite[Theorem 4.3]{sag} for details), we get that if $k<k'$ then there are $2^{k-1}$ partitions of size $k$ avoiding a $(3,k')$-pair but only $2^{k-1}-1$ partitions of size $k$ avoiding a $(3,k)$-pair.
Hence, no $(3,k)$-pair can be Wilf-equivalent to a $(3,k')$-pair where $k \neq k'$.

A detailed discussion of the above facts can be found in the paper of Jel\'{\i}nek, Mansour, and Shattuck \cite{jel}.
In fact, they also showed that all except one Wilf-equivalence class of $(3,k)$-pairs correspond to dominating equivalence classes of $C_k$.

\begin{enumerate}
\setlength{\itemsep}{0.2cm}
    \item The Wilf-equivalence class corresponding to the dominating equivalence class $E$ containing $(1,\ldots,1)$ consists of the following pairs:
    \begin{enumerate}
    \setlength{\itemsep}{0.15cm}
        \item $\{121,\tau_{121}(a)\}$ where $a \in E$,
        \item $\{112,\tau_{112}(a)\}$ where $a \in E$,
        \item $\{122,\tau\}$ where $\tau \in P_k(\{122\})$, and
        \item $\{123,\tau'\}$ where $\tau' \neq 1^k$.
    \end{enumerate}
    \item The Wilf-equivalence class corresponding to a dominating equivalence class $E$ not containing $(1,\ldots,1)$ consists of the following pairs:
    \begin{enumerate}
    \setlength{\itemsep}{0.15cm}
        \item $\{121,\tau_{121}(a)\}$ where $a \in E$, and
        \item $\{112,\tau_{112}(a)\}$ where $a \in E$.
    \end{enumerate}
    \item The pair $\{123,1^k\}$ is not Wilf-equivalent to any other $(3,k)$-pair.
\end{enumerate}

Hence, if we denote the number of dominating equivalence classes in $C_k$ as $\xi_k$, we get that $(3,k)$-pairs split up into $1+\xi_k$ Wilf-equivalence classes.
Hence, we now shift our focus to the dominating equivalence among compositions.

\section{Dominating equivalence}\label{sec3}

We first recall a few results from the paper of
Jel\'{\i}nek, Mansour, and Shattuck \cite{jel}.

\begin{lemma}[{\cite[Lemma 2.3]{jel}}]\label{lemshuf}
If $a=(a_1,\ldots,a_m)$ and $a'=(a_{\sigma(1)},\ldots,a_{\sigma(m)})$ where $\sigma$ is a permutation of $[m]$, then $a \overset{d}{\sim} a'$.
\end{lemma}

\begin{lemma}[{\cite[Lemma 2.4]{jel}}]\label{lemcomb11}
If $a=(a_1,\ldots,a_{m-1},2)$ and $a'=(a_1,\ldots,a_{m-1},1,1)$, then $a \overset{d}{\sim} a'$.
\end{lemma}

These lemmas imply that to each composition we can associate a unique $2$-free integer partition which is dominating equivalent to it.
This is done by first rearranging the terms in decreasing order and then replacing each $2$ in the composition by two $1$s.

\begin{example}
The composition $(3,2,5,1,1,3)$ is converted to the associated unique $2$-free integer partition as follows:
\begin{equation*}
    (3,2,5,1,1,3) \overset{d}{\sim} (5,3,3,2,1,1) \overset{d}{\sim} (5,3,3,1,1,1,1).
\end{equation*}
\end{example}

Hence, any two compositions that correspond to the same $2$-free integer partition are dominating equivalent.
We will now show that the converse is true as well, which was conjectured by Jel\'{\i}nek, Mansour, and Shattuck \cite{jel}.

\begin{theorem}\label{thm2free}
Two compositions $a$ and $a'$ are dominating equivalent if and only if they correspond to the same $2$-free integer partition.
\end{theorem}

The above theorem implies that $\xi_k$, i.e., the number of dominating equivalence classes in $C_k$, is $p(k)-p(k-2)$, where for any $n \geq 1$, $p(n)$ is the number of integer partitions of $n$.
The sequence $(\xi_k)_{k\geq 0}$ is listed in the OEIS \cite{oeis} as \href{http://oeis.org/A027336}{A027336}.
Hence the number of Wilf-equivalences classes of $(3,k)$-pairs is $p(k)-p(k-2)+1$.

To prove the theorem, we first need a few lemmas and definitions.

\begin{lemma}\label{lem1}
If $a \overset{d}{\sim} a'$ for two compositions $a$ and $a'$ of $m$ and $n$ respectively, then $m=n$.
\end{lemma}
\begin{proof}
Suppose to the contrary that $m$ and $n$ are distinct.
Without loss of generality, we assume $m<n$.
Then $D_m(a)=\{a\}$ whereas $D_m(a')$ is an empty set.
This contradicts the fact that $a \overset{d}{\sim} a'$ and hence proves our lemma.
\end{proof}

\begin{definition}
The composition $b=(b_1,\ldots,b_k)$ is said to \textit{h-dominate} the composition $a$ if $b$ dominates $a$ and either $k=1$ or $(b_1,\ldots,b_{k-1})$ does not dominate $a$.
For any positive integer $n$, $D^h_n(a)$ is the set of compositions of $n$ that h-dominate $a$.
We say $a$ and $a'$ are \textit{h-dominating equivalent}, written as $a \overset{h}{\sim} a'$, if for all positive integers $n$, $|D_n^h(a)|=|D_n^h(a')|$.
\end{definition}

We will now show that dominating equivalence is the same as h-dominating equivalence.
Before doing so, we set up some notations.
Let $c=(c_1,\ldots,c_k)$ and $c'=(c'_1,\ldots,c'_{k'})$ be compositions.
We will denote the composition $(c_1,\ldots,c_k,c'_1,\ldots,c'_{k'})$ by $(c,c')$.
Similarly, for any positive integer $s$, denote the composition $(c_1,\ldots,c_k,s)$ by $(c,s)$.

\begin{lemma}\label{lem2}
For any two compositions $a$ and $a'$, $a \overset{d}{\sim} a'$ if and only if $a \overset{h}{\sim} a'$.
\end{lemma}
\begin{proof}
Note that for any composition $a$,
\begin{equation}\label{eq1}
    D_n(a) = D_n^h(a)\ \sqcup\ \bigsqcup_{i=1}^{n-1} \{(c,c') \mid c \in D_i^h(a),\ c' \in C_{n-i}\}.
\end{equation}
This is obtained by finding, for each $b=(b_1,\ldots,b_k) \in D_n(a)$, the least $j \in [k]$ such that $(b_1,\ldots,b_j)$ dominates $a$.
This element $b$ of $D_n(a)$ belongs to exactly one of the $n$ sets on the right side.
If $j=k$, then $b \in D_n^h(a)$ or else $b \in \{(c,c') \mid c \in D_i^h(a),\ c' \in C_{n-i}\}$ where $i=b_1+\cdots+b_j$.

Suppose $a \overset{h}{\sim} a'$.
The equality \eqref{eq1} implies that for all positive integers $n$,
\begin{align*}
    |D_n(a)| &= |D_n^h(a)|\ +\ \sum_{i=1}^{n-1} |D_i^h(a)|\cdot|C_{n-i}|\\
    &= |D_n^h(a')|\ +\ \sum_{i=1}^{n-1} |D_i^h(a')|\cdot|C_{n-i}|\\
    &= |D_n(a')|.
\end{align*}
Hence we get $a \overset{d}{\sim} a'$.

Conversely, suppose $a \overset{d}{\sim} a'$.
We will prove $|D_n^h(a)|=|D_n^h(a')|$ for all $n \geq 1$ by induction on $n$.
We have $|D_1^h(a)|=|D_1(a)|=|D_1(a')|=|D^h_1(a')|$.
Let $n \geq 2$ and $|D_k^h(a)| = |D_k^h(a')|$ for all $k<n$.
Using \eqref{eq1}, we get
\begin{align*}
    |D_n^h(a)| &= |D_n(a)|\ -\ \sum_{i=1}^{n-1} |D_i^h(a)|\cdot|C_{n-i}|\\
    &= |D_n(a')|\ -\ \sum_{i=1}^{n-1} |D_i^h(a')|\cdot|C_{n-i}|\\
    &= |D_n^h(a')|.
\end{align*}
Hence we get $a \overset{h}{\sim} a'$, which completes the proof of the lemma.
\end{proof}

\begin{lemma}\label{lem3}
If $a$ and $a'$ are two compositions and $s$ is a positive integer such that $(a,s) \overset{h}{\sim} (a',s)$, then we have $a \overset{h}{\sim} a'$.
\end{lemma}
\begin{proof}
Let $a$ be a composition and $s$ be a positive integer.
Let $C_{k,s}$ denote those compositions of $k$ of the form $(c_1,\ldots,c_m)$ where $c_i<s$ for all $i \in [m-1]$ and $c_m \geq s$.
Note that
\begin{equation}\label{eq2}
    D^h_n((a,s)) = \bigsqcup_{i=1}^{n-1}\{(c,c') \mid c \in D^h_i(a),\ c' \in C_{n-i,s}\}.
\end{equation}
This is obtained by finding, for each $b=(b_1,\ldots,b_m) \in D^h_n((a,s))$, the least $j \in [m]$ such that $(b_1,\ldots,b_j)$ dominates $a$.
Since $b \in D^h_n((a,s))$, this would mean that $b_{j+1},\ldots,b_{m-1}$ are all strictly less than $s$ and $b_m \geq s$.
Note that $|C_{i,s}|=0$ for all $i < s$ and that $|C_{s,s}|=1$.
Hence an implication of the equality \eqref{eq2} is that for any positive integer $n$,
\begin{equation}\label{eq3}
    |D_{n+s}^h((a,s))| = |D_n^h(a)| + \sum_{i=s+1}^{n+s-1}|D_{n+s-i}^h(a)| \cdot |C_{i,s}|.
\end{equation}
Let $a'$ be a composition such that $(a,s) \overset{h}{\sim} (a',s)$.
We will prove $|D_n^h(a)|=|D_n^h(a')|$ for all $n \geq 1$ by induction on $n$.
If $|D_1^h(a)| \neq |D_1^h(a')|$, we would have exactly one of $a$ or $a'$ being the composition of $1$.
This would mean $(a,s)$ and $(a',s)$ are compositions of different numbers, which contradicts \cref{lem1} (since \cref{lem2} implies that $(a,s) \overset{d}{\sim} (a',s)$).
Let $n \geq 2$ and $|D_k^h(a)| = |D_k^h(a')|$ for all $k<n$.
From \eqref{eq3}, we get
\begin{align*}
    |D_n^h(a)| &= |D_{n+s}^h((a,s))| - \sum_{i=s+1}^{n+s-1}|D_{n+s-i}^h(a)| \cdot |C_{i,s}|\\
    &= |D_{n+s}^h((a',s))| - \sum_{i=s+1}^{n+s-1}|D_{n+s-i}^h(a')| \cdot |C_{i,s}|\\
    &= |D_n^h(a')|.
\end{align*}
Hence we get $a \overset{h}{\sim} a'$, which completes the proof of the lemma.
\end{proof}

\begin{corollary}\label{cor}
If $a$ and $a'$ are two compositions and $s$ is a positive integer such that $(a,s) \overset{d}{\sim} (a',s)$, then we have $a \overset{d}{\sim} a'$.
\end{corollary}
\begin{proof}
Combine \cref{lem2} and \cref{lem3}.
\end{proof}

We can now prove our main theorem.

\begin{proof}[Proof of \cref{thm2free}]
We have to show that if $a$ and $b$ are two compositions such that $a \overset{d}{\sim} b$, then they correspond to the same $2$-free integer partition.
By \cref{lem1}, such $a$ and $b$ are compositions of the same number $n$.
We will prove the theorem by induction on $n$, the case $n=1$ being trivial.
Let $n \geq 2$ and suppose the statement is true for all numbers less than $n$.
Suppose $a \overset{d}{\sim} b$ are two compositions of $n$.
Using \cref{lemshuf} and \cref{lemcomb11}, we can assume that $a_1 \geq \cdots \geq a_m$ and $a_{m-1} \geq 2$ and similarly that $b_1 \geq \cdots \geq b_k$ and $b_{k-1} \geq 2$ (this can be done by reordering the terms in decreasing order and then changing any pair of $1$'s to a $2$).

We will now compute $|D_{n+1}(a)|$.
If $a_m \geq 2$, we get $|D_{n+1}(a)|=2m+1$.
This is because any composition in $D_{n+1}(a)$ is obtained by either
\begin{enumerate}
    \item inserting the term $1$ before $a_1$, between $a_i$ and $a_{i+1}$ for some $i \in [m-1]$, or after $a_m$, or
    \item adding $1$ to some term of $a$.
\end{enumerate}
Since all terms of $a$ are greater than $1$ in this case, each method of obtaining a composition of $D_{n+1}(a)$ described above results in a different composition.

When $a_m=1$, using the same logic as above but noting that adding $1$ before $a_m$ and adding $1$ after $a_m$ both result in the same composition, we get $|D_{n+1}(a)|=2m$.

Using similar arguments for $b$ and since we must have $|D_{n+1}(a)|=|D_{n+1}(b)|$, we get $k=m$, i.e., $a$ and $b$ have the same number of terms.

Note that if we show $a_m=b_m$, \cref{lem3} and the induction hypothesis would imply that $a$ and $b$ have the same corresponding $2$-free integer partition.
On the contrary, let $a_m \neq b_m$.
Without loss of generality we can assume that $a_m < b_m$.
We will show that this implies that $|D_{n+a_m}(a)| < |D_{n+a_m}(b)|$, which is a contradiction to $a \overset{d}{\sim} b$.

A \textit{constructive pair} is a pair $(c,P)$ where $c=(c_1,\ldots,c_l)$ is a composition of $a_m$ and $P$ is a sequence that consists of the numbers $1,2,\ldots,l$ in order and $m$ boxes such that each box contains at most one of the numbers.
To each such pair $(c,P)$ we associate the composition in $D_{n+a_m}(b)$ obtained by replacing each unboxed number $i$ in the sequence $P$ by $c_i$ and the $j^{th}$ box by $b_j$ if it is empty and by $b_j+c_i$ if it contains the number $i$.
We will call this composition $b(c,P)$.

\begin{example}\label{egconst}
Suppose $b=(8,6,6,5)$, $a_m=4$, $c=(1,2,1)$, and $P=1\ \framebox{\color{white}2}\ \framebox{2}\ 3\ \framebox{\color{white}2}\ \framebox{\color{white}2}$.
Then the composition $b(c,P)$ is $(1,8,8,1,6,5)$.

The construction of $b(c,P)$ from $b$ be thought of visually as follows:
Think of any composition as towers of boxes.
The boxes in the sequence $P$ can be thought of as the top view of $b$ and the number $i$ represents where to add $c_i$ boxes to $b$.
The construction in the example is shown in \Cref{fig:my_label}.
\end{example}

\begin{figure}[H]
    \centering
    \begin{tikzpicture}[scale=0.5]
        \draw (0,0) rectangle (0+1,0+1);
        \draw (0,1) rectangle (0+1,1+1);
        \draw (0,2) rectangle (0+1,2+1);
        \draw (0,3) rectangle (0+1,3+1);
        \draw (0,4) rectangle (0+1,4+1);
        \draw (0,5) rectangle (0+1,5+1);
        \draw (0,6) rectangle (0+1,6+1);
        \draw (0,7) rectangle (0+1,7+1);

        \draw (1,0) rectangle (1+1,0+1);
        \draw (1,1) rectangle (1+1,1+1);
        \draw (1,2) rectangle (1+1,2+1);
        \draw (1,3) rectangle (1+1,3+1);
        \draw (1,4) rectangle (1+1,4+1);
        \draw (1,5) rectangle (1+1,5+1);

        \draw (2,0) rectangle (2+1,0+1);
        \draw (2,1) rectangle (2+1,1+1);
        \draw (2,2) rectangle (2+1,2+1);
        \draw (2,3) rectangle (2+1,3+1);
        \draw (2,4) rectangle (2+1,4+1);
        \draw (2,5) rectangle (2+1,5+1);

        \draw (3,0) rectangle (3+1,0+1);
        \draw (3,1) rectangle (3+1,1+1);
        \draw (3,2) rectangle (3+1,2+1);
        \draw (3,3) rectangle (3+1,3+1);
        \draw (3,4) rectangle (3+1,4+1);
        \node at (6,3) {$\xrightarrow{(c,P)}$};
        \node at (2,-1) {$\textcolor{white}{(P}b\textcolor{white}{P)}$};
    \end{tikzpicture}
    \begin{tikzpicture}[scale=0.5]
        \draw[blue] (-1,0) rectangle (-1+1,0+1);
        \node at (-0.5,0.5) {\color{blue}$c_1$};
        \draw (0,0) rectangle (0+1,0+1);
        \draw (0,1) rectangle (0+1,1+1);
        \draw (0,2) rectangle (0+1,2+1);
        \draw (0,3) rectangle (0+1,3+1);
        \draw (0,4) rectangle (0+1,4+1);
        \draw (0,5) rectangle (0+1,5+1);
        \draw (0,6) rectangle (0+1,6+1);
        \draw (0,7) rectangle (0+1,7+1);

        \draw (1,0) rectangle (1+1,0+1);
        \draw (1,1) rectangle (1+1,1+1);
        \draw (1,2) rectangle (1+1,2+1);
        \draw (1,3) rectangle (1+1,3+1);
        \draw (1,4) rectangle (1+1,4+1);
        \draw (1,5) rectangle (1+1,5+1);
        \draw[blue] (1,6) rectangle (1+1,6+2);
        \node at (1.5,7) {\color{blue}$c_2$};

        \draw[blue] (2,0) rectangle (2+1,0+1);
        \node at (2.5,0.5) {\color{blue}$c_3$};
        \draw (2+1,0) rectangle (2+1+1,0+1);
        \draw (2+1,1) rectangle (2+1+1,1+1);
        \draw (2+1,2) rectangle (2+1+1,2+1);
        \draw (2+1,3) rectangle (2+1+1,3+1);
        \draw (2+1,4) rectangle (2+1+1,4+1);
        \draw (2+1,5) rectangle (2+1+1,5+1);

        \draw (3+1,0) rectangle (3+1+1,0+1);
        \draw (3+1,1) rectangle (3+1+1,1+1);
        \draw (3+1,2) rectangle (3+1+1,2+1);
        \draw (3+1,3) rectangle (3+1+1,3+1);
        \draw (3+1,4) rectangle (3+1+1,4+1);
        \node at (3,-1) {$b(c,P)$};
        \end{tikzpicture}
    \caption{Construction corresponding to the constructive pair given in \cref{egconst}}
    \label{fig:my_label}
\end{figure}

Using any occurrence of $b$ in a composition $b'$ in $D_{n+a_m}(b)$, we can obtain a constructive pair $(c,P)$ such that $b(c,P)=b'$.
This is done as follows:
Suppose $b'=(b'_1,\ldots,b'_p)$ is a composition in $D_{n+a_m}(b)$ and $1 \leq i(1) < i(2) < \cdots < i(m) \leq p$ is an occurrence of $b$.
Construct $P$ from $b'$ by first replacing $b'_j$ by a ball if $j \notin \{i(1),\ldots,i(m)\}$.
Then for any $j \in [m]$, replace $b'_{i(j)}$ by an empty box if $b'_{i(j)}=b_j$ and by a box containing a ball if $b'_{i(j)}>b_j$.
Suppose there are $l$ balls in $P$, we replace them by the numbers $1,\ldots,l$ in order.
Now for each $i \in [l]$,
\begin{enumerate}
    \item set $c_i$ as $b'_j$ if $i$ is an unboxed number and the $j^{th}$ term in $P$, and
    \item set $c_i$ as $b'_{i(j)}-b_j$ if $i$ is a boxed number and the $i(j)^{th}$ term of $P$.
\end{enumerate}
The above procedure might be more clear when compositions are thought of as towers of boxes.
For example, in \Cref{fig:my_label}, if the black boxes in $b(c,P)$ are used as the occurrence of $b$ in the above procedure, we obtain the constructive pair $(c,P)$ of \cref{egconst}.

It is clear that if $b'=b(c,P)$ for some constructive pair $(c,P)$, then there is a canonical occurrence of $b$ in $b'$ formed using the original copy of $b$ to which terms were added in the construction.
Using this canonical occurrence in the above procedure we can get back the original constructive pair $(c,P)$.
However, since $b$ has all terms strictly greater than $a_m$, there is a unique occurrence of $b$ in any composition $b'$ in $D_{n+a_m}(b)$.
Hence, there is a unique constructive pair $(c,P)$ such that $b(c,P)=b'$ for any composition $b'$ of $D_{n+a_m}(b)$.
Therefore, the set $D_{n+a_m}(b)$ is in bijection with the set of constructive pairs.

Since $a$ also has $m$ terms, the same method can be used to associate a composition in $D_{n+a_m}(a)$ to a constructive pair $(c,P)$, and we call this composition $a(c,P)$.
Just as for $b$, any composition $a' \in D_{n+a_m}(a)$ can be obtained from $a$ using a constructive pair.
However, $a((a_m),P_1)=a((a_m),P_2)$ where $P_1$ has $1$ between the $(m-1)^{th}$ and $m^{th}$ box and $P_2$ has $1$ after the $m^{th}$ box.
Hence we get $|D_{n+a_m}(a)| < |D_{n+a_m}(b)|$, which completes the proof.
\end{proof}

\bibliographystyle{abbrv}