\section{Conclusion}
This report describes our systems for a shared task of evaluating the rationales of amateur investors at FinNLP-2022. From the experimental results in pairwise comparison, we conclude that few-shot prompted instruction-based language models can work reasonably well in low resource settings with minimal training efforts but might need quite accurate data from sources since using translated data seems not to perform well. From the experimental results in unsupervised ranking, financial language models perform well and Bayesian-fitting the lexicons helps improve the performance. Also, the heuristic function design needs to differ between MPP and ML.
\section*{Limitations}
We only sampled a relatively small portion of models and draw conclusions. We also conducted experiments only on one dataset for evaluating the rationales of amateur investors. Besides, the dataset is in Chinese with English translation using Google Translate. Lots of our methods rely on the translated data.

Because we are limited to only three submissions, we don't know how each feature set contributes to the score. There were no ablations. However, the shared task organizers released the test data with ground truths in private.

The authors are self-affiliated and do not represent any entities. The authors also participated in the shared task under many severe unattended local personal criminal events in their home countries. There might be some unintentional errors and physical limitations based on these unlawful interruptions. Even at the time of drafting this report, the authors suffer from unknown toxin flumes spraying into their places. We want to participate in the shared task because it is fun and educational. We apologize for any errors in this report. We tried our best.

\section*{Ethics Statement}
Scientific work published at EMNLP 2022 must comply with the \href{https://www.aclweb.org/portal/content/acl-code-ethics}{ACL Ethics Policy}. We, the authors, hope the intended uses of our systems are for peace, well-being, and social good only. No harm.

\section*{Acknowledgments}
We would like to thank anonymous reviewers for their constructive feedback and additional experiment suggestions.

\section{Experimental Results}
In our experiments, most of our submissions (except T5-small) are intuition-based heuristics, and we did not even measure neither their training nor validation performance at all during the competition. We did not use any data augmentation techniques.
\begin{table}[t]
\centering
\caption{MPP and ML accuracies of our models in pairwise comparison test data. (The numbers in subscript are submission rankings on the leaderboard. The symbol $\dagger$ denotes a top-3 performance.)}
\vspace{-0.5\baselineskip}
\label{st1table}
\begin{tabular}{|l|c|c|}
\hline
 Model    &  MPP acc. & ML acc.\\
\hline
T5-small  & $47.13_{16}$ & $45.98_{10}$ \\
InstructGPT-zh  & $\mathbf{48.28_{14}}$ & $\mathbf{54.02_3}\dagger$ \\
InstructGPT-en  & $47.13_{16}$ & $41.38_{13}$\\
\hline
FinNLP-22 best & $62.07$& $59.77$ \\
\hline
\end{tabular}
\end{table}

\begin{table}[t]
\centering
\caption{Average MPP and ML from top $10\%$ posts of our models in unsupervised ranking test data. (The numbers in subscript are submission rankings on the leaderboard. The symbol $\dagger$ denotes a top-3 performance and the symbol $\ddagger$ denotes the score beats the baseline.)}
\vspace{-0.5\baselineskip}
\label{st2table}
\begin{tabular}{|l|c|c|}
\hline
 Model    &  avg. MPP & avg. ML \\
 \hline
Stylistic baseline & $17.61\%$& $-2.46\%$ \\
\hline
base-1  & $22.53\%_{3}\dagger\ddagger$ & $\mathbf{-7.80\%_{11}}$ \\
bayesdcm-2  & $\mathbf{24.39\%_{1}\dagger\ddagger}$ & $-13.04\%_{16}$ \\
multinomial-3  & $23.76\%_{2}\dagger\ddagger$ & $-12.33\%_{15}$\\
\hline
FinNLP-22 best & $24.39\%$ (ours)& $-2.46\%$ \\
\hline
\end{tabular}
\vspace{-1\baselineskip}
\end{table}

\subsection{Pairwise Comparison}
The experimental results in Table.\ref{st1table} show that the OpenAI InstructGPT language model few-shot trained on Chinese data works best in our submissions, ranking 3rd on the maximal loss (ML) pairwise accuracy, even better than instead training on the Google translated English data by a large margin, where the English few-shot trained InstructGPT model even performs worse than an instruction-based T5-small model finetuned on the English data. However, all instruction-based submissions do not perform well on the maximal potential profit (MPP) pairwise accuracy where there are more data and learning signals, nonetheless, the Chinese few-shot trained InstructGPT model still performs best in our setting.

We additionally split the training data into a held-out train/val split and evaluated our methods on the val split in Table \ref{st3table}. The results are a bit different since the English version of the InstructGPT works better. However, we did not hope for an accurate cross-validation estimation given a small amount of data. Using leave-one-out validation (LOOCV) or $k$-fold cross-validation with a high value of $k$ can produce a better estimation but they are costly. We might be able to generate more data pairs, but we decided to keep the same setting.

\begin{table}[t]
\centering
\caption{Additional experiments on using our pairwise comparison methods on a held-out train/val split (ratio=$0.3$). The evaluation metric is accuracy.}
\vspace{-0.5\baselineskip}
\label{st3table}
\begin{tabular}{|l|c|c|}
\hline
 Model    &   MPP acc. &  ML acc.\\
\hline
T5-small  & 0.4833 & \textbf{0.6000}\\
InstructGPT-zh  & 0.4667 & 0.4167\\
InstructGPT-en  & \textbf{0.6167} & 0.4667 \\
\hline
\end{tabular}
\vspace{-1\baselineskip}
\end{table}

\subsection{Unsupervised Ranking}
For the unsupervised ranking task, we utilized many language models, including many financial-specific ones, and Bayesian lexicons unsupervisely learned on both Chinese and English words. All of our submissions rank best in the MPP ranking, from 1st to 3rd in this task. However, they all do not perform well for the ML scoring. Therefore, both MPP and ML scores need different treatments substantially since we treated MPP and ML using the same formula. Our only difference is the treatment of market sentiment lexicons. We feel that the sentiment features, or mostly semantic features, might be negatively correlated, weakly correlated, or even uncorrelated with ML because the stylistic baseline performs best, and our base submission performs better than our Bayesian lexicon submissions.

We conducted additional experiments on unsupervised ranking by using the whole training set of the pairwise comparison data. We compared all posts using our scoring functions in Table \Ref{st4table}. The results show not much difference among our methods. When we tried to evaluate using the pairwise comparison accuracy, the results show no difference ($0.545$ MPP comparison acc. and $0.525$ ML comparison acc.) as our methods were not designed for that.

\section{Introduction}
Evaluating the rationals of amateur investors (ERAI) \cite{chen2021evaluating, chen2021opinion} is a shared task on evaluating social media opinions on the topic of investments and whether they are going to be useful or not. Mining high-quality opinions by inspecting their supporting rationales might utilize the wisdom of the crowd on social media. Previous work \cite{chen2021evaluating} proposes stylistic and semantic features to filter out noisy crowd opinions which may not be high-quality and profitable. There are two settings in this ERAI shared task, pairwise comparison and unsupervised ranking. These settings sort out the opinions based on two metrics, higher maximal potential profit (MPP) and lower maximal loss (ML). In pairwise comparison, two posts are given with a binary label whether the MPP and ML of the first post are more or less than the second post. In unsupervised ranking, the goal is to filter and keep the top 10\% posts based on MPP and ML given a set of unranked posts.

For pairwise comparison, our best submission ranks $3^{rd}$ on the maximal loss (ML) pairwise accuracy on the leaderboard\footnote{\url{https://sites.google.com/nlg.csie.ntu.edu.tw/finnlp-2022-emnlp/erai-shared-task}}. For unsupervised ranking, our best submission ranks $1^{st}$ on the maximal potential profit (MPP) ranking. The codes for our systems are open-sourced and available at our GitHub repository\footnote{\url{https://github.com/perathambkk/finnlp_erai_shared_task_emnlp2022}}.
\section{Models}
\subsection{Pairwise Comparison}
For pairwise comparison, we utilized instruction-based models based on T5-small \cite{raffel2020exploring} and OpenAI InstructGPT language models \cite{ouyang2022training} in a few-shot prompt-based setting \cite{brown2020language}.

\subsubsection{T5}
T5 is an encoder-decoder language model which was trained by treating every text processing problem as a ``text-to-text" problem to unify NLP tasks using only a single model, loss function, hyperparameter set, etc. The input texts will be encoded and the T5 decoder will decode them. Specifically, T5 was unsupervised-pretrained by denoising masked inputs on the ``Colossal Clean Crawled Corpus" (C4) dataset, Common Crawled from web scraping. Then, T5 can be further supervised finetuned using the ``text-to-text'' format and T5's decoder will be in the teacher forcing mode where the decoder will be trained using input and a right-shifted target sequence. T5 architecture is pretty much the same as the vanilla transformer \cite{vaswani2017attention} except removing the LayerNorm bias term, placing the LayerNorm outside the residual path, and using a relative position embedding \cite{shaw-etal-2018-self}.

In the T5 paper, the authors state that T5 can be specified which task it should perform by adding a task-specific textual prefix to the original input sequence before feeding into the model. Therefore, we take T5-small as an instruction-based model using the input prompt template, `post1 : \%s post2 : \%s  </s>', where the \%s contains texts from the corresponding post, and the output prompt, `maximal potential profit (MPP) : \%s maximal loss (ML) : \%s </s>', where the \%s contains the MPP and the ML corresponding labels accordingly. This is similar to the baseline system in the FLUTE figurative language understanding dataset paper \cite{chakrabarty2022flute}, however, in our case, the T5-small is expected to jointly predict both MPP and ML in one forward pass of preparing the probability tensor. We use top-p sampling for text generation \cite{holtzman2019curious}.

\subsubsection{OpenAI's InstructGPT}
The OpenAI API has many variants of InstructGPT language models based on the GPT-3 autoregressive language model to conveniently perform various NLP tasks with the prompt library. The InstructGPT was trained with a human-in-the-loop style and is claimed to be better at following instructions, more truthful, and less toxic than the GPT-3. In this shared task, we engineered the prompts for InstructGPT using a few-shot learning setting, as in the GPT-3 paper \cite{NEURIPS2020_1457c0d6}, where few data instances were given from the target task/domain. Each data instance will become a prompt as `post1 : d[`post1'] post2: d[`post2'] > maximal potential profit (MPP)| \%s\# maximal loss (ML)| \%s.', where d[`post1'] and d[`post2'] are texts from the corresponding post. Then, we append the query we want to predict MPP and ML as just a truncated template, `post1 : d[`post1'] post2: d[`post2'] >', and let the language model generate the rest.

We use the `text-davinci-002' model and randomly construct those few-shot prompts where each prompt will be a length of around $4,000$ because of the API token length limit. We use the same setting and the model pipeline for both of our submissions $2$ and $3$ where we use the Chinese posts as d[`post1'] and d[`post2'] for our submission $2$ and the Google-translated English posts for our submission $3$. By this we mean, for example, we use the same tokenizer for Chinese and English. Therefore, these systems are very simple and to-go prompt-based systems. We had done very minimal parameter tuning to the model, only prompt engineering. For a survey in prompt-based systems, please consider \cite{liu2021pre}.

\subsection{Unsupervised Ranking}
For unsupervised ranking, we utilized many financial and general language models and Bayesian lexicons in both Chinese and English.

\subsubsection{Base Model}
Our first submission, our base model, consists of a stylistic length feature \cite{zong-etal-2020-measuring} derived from the opinion (sub)word lengths segmented using the `hfl/chinese-bert-wwm-ext' tokenizer \cite{cui2021pre}, prediction scores from FinBERT-FLS \cite{huang2020finbert}, a professional lexicon count from FinProLex \cite{chen2021evaluating}, and a market sentiment lexicon count from NTUSD-Fin \cite{chen2018ntusd}.

In the measuring forecasting skill from text paper \cite{zong-etal-2020-measuring}, the authors observe various linguistics phenomena indicating that skilled forecasters tend to write significantly longer justifications because of more rationale. For example, skilled forecasters also provide less readability, because of the usage of more complex languages, and less emotion, because of the usage of less emotional languages as neutral sentiments. Moreover, skilled forecasters tend to use more cardinal numbers, prepositions, and nouns. They tend to use fewer verbs and pronouns. Therefore, in this base model, we just stick with the lengths of justifications as our simplest skill indicator.

FinBERT \cite{huang2020finbert} is essentially BERT \cite{devlin-etal-2019-bert} customized for financial texts, pretrained on corporate filings, analyst reports, and earnings conference call transcripts, which differ from normal texts in both vocabulary and writing style. In the FinBERT paper, FinBERT outperforms all other methods, Loughran McDonald lexicon, and machine learning algorithms, especially in negative financial sentiment prediction, when finetuned for financial sentiment analysis (FinBERT-tone). FinBERT was finetuned in two additional tasks, labeling environment, social, and governance (ESG) discussions and labeling forward-looking statements (FLS), from firms’ corporate social responsibility (CSR) reports and management discussion and analyses (MD\&As) textual sentences. For this base model, we sum and normalize the prediction logit outputs from FinBERT-FLS classes, $\{FLS, NON\_FLS, NOT\_FLS\}$ on each sentence of the textual inputs as our scores.

FinProLex \cite{chen2021evaluating} is a Chinese financial lexicon derived from Bloomberg Terminal and PTT Stock Taiwanese social media platform, containing $5,162$ tokens from professional analysts' reports and social media posts paired with expertise scores. FinProLex uses Point-wise Mutual Information (PMI), as in \cite{turney-2002-thumbs, li-shah-2017-learning}, to measure the association strengths between a word and either the positive or negative lexicon. The formula of the expert-like score (ELScore) of a given word $w$ is as follows:
\begin{multline}
    ELScore_w = PMI(w, analyst) \\
    - PMI(w, amateur),
\end{multline}
\begin{multline}
   ELScore_w =\log_{2}\frac{p(w, analyst)}{p(w)p(analyst)} \\- \log_{2}\frac{p(w, amateur)}{p(w)p(amateur)},
\end{multline}
where $analyst$ and $amateur$ are labels of whether a given word is from an analyst report or an amateur post. This is the difference value between the PMI scores measuring how much a term is associated with either analyst or amateur documents. This is similar to the term's sentiment score ($S_{PMI}$) \cite{li-shah-2017-learning} which is
\begin{equation}\label{eq:spmi}
S_{PMI} = PMI(w, bullish) - PMI(w, bearish).
\end{equation}

FinProLex tends to include hard words, complex semantics, noun phrase modifiers, content words, transition words, personal pronouns, and negative words as experts tend to use most of them, except personal pronouns and negative words which are used more by amateurs, based on the paper findings that can be summarized as experts tend to evaluate pricing and valuations while amateurs tend to predict the stock movements.

NTUSD-Fin is an English lexicon for market sentiment analysis from StockTwits, containing $8,331$ words, $112$ hashtags, and $115$ emojis. We used their market sentiment scores which are also computed essentially from equation (\ref{eq:spmi}).

We aggregated the scores to predict MPP and ML using these heuristic functions ($\operatorname{base-1}$),
\begin{multline}
MPP = len + FinProLex + |FinWord>0| \\+ (FLS + 0.5  \operatorname{NON\_FLS} - \operatorname{NOT\_FLS}),
\end{multline}
\begin{multline}
ML = len + FinProLex + |FinWord<0| \\+ (FLS + 0.5 \operatorname{NON\_FLS} - \operatorname{NOT\_FLS}).
\end{multline}

We simply used a weighted sum as our heuristic function. We grouped similar scores together. $len + FinProLex$ are stylistic features where we put an equal weight of $1$ for each of them. We used $FinWord$ as a switch feature for either MPP or ML that would behave differently because of the market sentiment based on our belief. MPP posts should be from a bullish market while ML posts should instead be from a bearish market. $FLS$ has $3$ different class scores so we weighted $1$ for a positive class, $0.5$ for a less positive one, and $-1$ for a negative class. The weights are just our rule-of-thumb (make-up numbers that we felt they made sense solely from our intuitions).

It is like trying to intuitively come up with a good feature weighting number for a Maximum Entropy (MaxEnt) model. From our heuristic functions, we just down-weighted some scores and specify some negative interactions. We did not normalize the weighting into probabilities but the ranking should be the same anyway. Most weightings are uniformly the same number.
\begin{itemize}
\item If a score should positively correlate with the target, we should give a high weight.
\item If a score should weakly correlate with the target, we should give a low weight.
\item If a score should negatively correlate with the target, we should give a high negative weight.
\item For the rest that we are not certain of, they should retain a maximum entropy (uniformity).
\item These heuristics can be estimated with intuitions and give an intuitive unsupervised aggregated scoring function.
\end{itemize}
 However, we submitted the same function for MPP and ML to get a sense of using the same strategy for both bullish and bearish markets.

\subsubsection{Bayesian Lexicons}
Next, we added Bayesian lexicons \cite{eisenstein2017unsupervised} (by fitting FinProLex and NTUSD-Fin), FinBERT-tone \cite{huang2020finbert}, (fitted) Loughran-McDonald financial sentiment lexicon (LM) \cite{loughran2011liability} and Part-of-Speech (POS) features \cite{zong-etal-2020-measuring} into the score aggregators. We would like to note that these lexicons are not multi-word (only unigrams) so they are not expected to be able to handle negations except the creators of those lexicons had made them handle some kind of negations, like in the LM lexicon Fin-Pos list. The authors use bigram to quadgrams counts when that bigram to quadgram follows some negation patterns. Our second and third submissions differ in the normalization of scores and Bayesian lexicon variants.

Loughran-McDonald financial sentiment lexicon (LM) \cite{loughran2011liability} was created because the Harvard Psychosociological Dictionary, specifically, the Harvard-IV-4 TagNeg (H4N) file, does not perform well in financial and accounting domains. Lots of Harvard dictionary negative words are not negative in finance.

Bayesian lexicon learns predictive weights for each word in a lexicon using a method-of-moments estimator from co-occurrence statistics without any labels as a special case of multinomial Na\"ive Bayes. For the second submission, we use the Dirichlet Compound Multinomial likelihood to reduce effective counts for repetitive words. For the third submission, we use the multinomial likelihood model. For example, when we fitted the LM lexicon using the pairwise comparison data, we gave $0.02626$ to `good', $0.00501$ to `optimistic', and $0.00278$ to `highest'. For LM negative words, we gave $0.00243$ to `decline', $0.00234$ to sharply, and $0.00186$ to `difficult'.

Our POS features are motivated by the measuring forecasting skill from text paper. We simply counted cardinal numbers, nouns, and verbs from Chinese jieba segmented texts. Then, these counts were normalized into the range of $[0, 1]$.

For these submissions, we sum and normalize the prediction logit outputs from FinBERT-tone classes, $\{pos\_tone, neg\_tone\}$ on each sentence of the textual inputs as our scores.

For the second submission, we aggregated the scores to predict MPP and ML using these heuristic functions ($\operatorname{bayesdcm-2}$),
\begin{multline}
MPP = len + FinProLex + |FinWord>0| \\ + (FLS + 0.5 \operatorname{NON\_FLS} - \operatorname{NOT\_FLS}) \\+ (pos\_tone - neg\_tone + LM) \\+ (nouns + cards - verbs),
\end{multline}
\begin{multline}
ML = len + FinProLex + |FinWord<0| \\+ (FLS + 0.5 \operatorname{NON\_FLS} - \operatorname{NOT\_FLS}) \\+ (pos\_tone - neg\_tone + LM) \\+ (nouns + cards - verbs).
\end{multline}

For the third submission, we aggregated the scores to predict MPP and ML using these heuristic functions ($\operatorname{multinomial-3}$),
\begin{multline}
MPP = 0.5(len + FinProLex) \\+ 0.33(FLS + 0.5 \operatorname{NON\_FLS} - \operatorname{NOT\_FLS})\\+ 0.33(pos\_tone - neg\_tone + LM) \\+ 0.33(nouns + cards - verbs)\\+ |FinWord>0|,
\end{multline}
\begin{multline}
ML = 0.5(len + FinProLex) \\+ 0.33(FLS + 0.5 \operatorname{NON\_FLS} - \operatorname{NOT\_FLS})\\+ 0.33(pos\_tone - neg\_tone + LM) \\+ 0.33(nouns + cards - verbs)\\+ |FinWord<0|.
\end{multline}
In these functions, we tried to group and reweigh the scores as normalization. If two or more scores mean the same thing, we might double count.


\section{Introduction}\vspace*{-1.5ex}
	In coherent optical communications, there is a great demand for flexible receiver algorithms which operate blindly and adapt to varying channel conditions. Operating receivers without known pilot symbols increases the effective data rate, which can be used to increase the forward error correction (FEC) overhead (and hence the robustness) or the net data throughput. However, the design of blind receivers is challenging, especially with modern higher-order modulation formats and \ac{PCS}.
	In particular, the bootstrapping of blind algorithms, especially adaptive equalizers, is a major hurdle which has not been widely studied in the literature and on which we focus in this work.

	A widely-used algorithm for bootstrapping blind equalization is the \ac{CMA}~\cite{godard}. Although designed for modulation formats with constant signal amplitude, it also converges for multi-amplitude formats such as $M$-QAM, where its criterion is ill-matched.
	There are variants of the \ac{CMA} for multi-amplitude signals, e.g., the \ac{RDE}~\cite{ready_RDE}, which, however, have worse convergence characteristics at startup, when the correct radii cannot be estimated reliably. Hence, the standard \ac{CMA} is often used for bootstrapping, and the algorithms are switched or the amount of radii is increased once the \ac{CMA} has converged~\cite{guiomar_fully_blind_2015,Savory_AS_2010digital}.
	However, the \ac{CMA} struggles with \ac{PCS} formats~\cite{zervas1991effects} or at certain working points. In this work, we compare different implementation variants of the \ac{CMA} at those critical points with the recently-proposed \ac{VAE} based equalizer~\cite{lauinger2022blind}.
	The latter explicitly considers the \textit{a~priori} density of the transmitted symbols, which can prevent the equalizer to converge to invalid constellations.

    \vspace*{-1ex}

\section{Blind Equalization} \vspace*{-1.5ex}
\begin{wrapfigure}[10]{R}{0.37\textwidth}
	\vspace*{-7ex}
	\raggedright
	\includegraphics{tikz_ext/System_VAELE_DP_new.pdf} \\%
	\includegraphics{tikz_ext/System_CMA_DP.pdf} \\%
	\vspace*{-2.5ex}
	\caption{Sketch of VAE (top) and CMA (bottom).} \label{Eq_structure}
\end{wrapfigure}
The widely-used equalizer for bootstrapping a coherent optical transmission system is the \ac{CMA}~\cite{Savory_AS_2010digital}. In its common implementation~\cite{godard}, the distance of the equalized signal amplitude to the kurtosis of the transmitted constellation is minimized by updating the equalizer taps \emph{symbol-wise}. By nature, the \ac{CMA} is phase-insensitive and requires an additional \ac{CPE}, where we apply the blind \emph{Viterbi-Viterbi} algorithm~\cite{viterbi1983nonlinear} with a sufficient averaging over 501~symbols. Furthermore, there are variants of the CMA applying batch-wise updating via \ac{GD}, e.g., as proposed in~\cite{crivelli2014}, which we denote as \emph{CMAbatch}.

In comparison, we analyze the recently proposed \ac{VAE} based equalizers~\cite{lauinger2022blind}, which use the \ac{ELBO} to approximate \ac{ML} channel estimation and, as a byproduct, equalize the received signal.
The structures of both equalizer types are shown in Fig.~\ref{Eq_structure}.
All algorithms update a similar linear equalizer block consisting of a complex-valued $2\times2$ butterfly structure with \ac{FIR} filters of length $M_{\mathrm{eq}}$,
and we use a maximum \textit{a~posteriori} based soft demapper as in \cite{cho2019probabilistic}.
In fact, the \ac{VAE} does not require a \ac{CPE}\footnote{Note that in the presence of laser phase noise, an additional \ac{CPE} might be necessary if the VAE cannot adapt rapidly enough.}
and jointly updates---along with the equalizer taps---a similar but separate butterfly structure $\bm{h}_\text{est}$ of $M_{\mathrm{est}}$-long \ac{FIR} filters for channel estimation. We use an Adam optimizer to update either batch-wise (\emph{VAEbatch}) or with the \emph{VAEflex} scheme as proposed in~\cite{lauinger2022blind}.
The latter is based on batch-wise processing, but instead of equalizing the whole batch per iteration, it only equalizes $N_\text{flex} <\Nb$~symbols. In the next iteration, it equalizes the next $N_\text{flex}$~symbols, so it re-uses $\Nb -N_\text{flex}$~symbols per iteration which boosts convergence. Further details can be found in~\cite{lauinger2022blind}. We also propose the \emph{CMAflex} scheme, which applies the same scheme based on the CMAbatch.

The batch-wise schemes consider a batch of $\Nb = 200$~symbols with learning rates of $\epsilon_{\mathrm{lr}}=\num{1.2e-4}$ (CMAbatch) and $\epsilon_{\mathrm{lr}}=\num{2e-3}$ (VAEbatch); the CMA uses $\epsilon_{\mathrm{lr}}=\num{8e-4}$. The CMAflex/VAEflex schemes equalize $N_\text{flex} = 10$~symbols while processing $\Nb = 100$~symbols per iteration, and use {$\epsilon_{\mathrm{lr}}=\num{4.5e-5}$}~(CMAflex) and $\epsilon_{\mathrm{lr}}=\num{2e-3}$ (VAEflex). All hyperparameters are optimized individually for each scheme to reduce the percentage of failed runs during bootstrapping.

\begin{figure} [!b]
		\centering
		\vspace*{-2ex}
		\includegraphics{tikz_ext/CD_full_v6.pdf}
		\vspace*{-4.5ex}
		\caption{Results for different residual \ac*{CD}: the top plots correspond to unshaped $64$-QAM ($\mathcal{H}=6$), the bottom plots to \ac*{PCS}-$64$-QAM ($\mathcal{H}=5.73$). The inset constellation diagrams show the equalized signal after the VAEbatch at $L_{\mathrm{cd}}=\SI{1.5}{\kilo\meter}$. The plots on the right show the mean frame index $\bar{k}$ at which a threshold of $\text{BMI}_\text{thr}=5$ and $\text{BMI}_\text{thr}=4.8$ (for the \ac*{PCS}, respectively) is reached.  }

		\label{Results:Lcd_full}
\end{figure}
\vspace*{-1ex}
\section{Simulation Environment}\vspace*{-1.5ex}
\def{\mathrm{j}}{{\mathrm{j}}}
To analyze the bootstrapping of the different equalizers, we adapt the simulation environment of \cite{lauinger2022blind} (see also references therein) and use the same averaging scheme.
Precisely, we transmit uniform $64$-QAM and \ac{PCS}-$64$-QAM with an entropy of $\mathcal{H}=5.73$ bits per symbol (and a Maxwell-Boltzmann distribution), assume a time-invariant channel during the bootstrapping phase, and model the fiber by the linear frequency domain channel matrix \vspace*{-0.7ex}
\begin{align*}
	\bm{H}\lr{f} &= \bm{R} \begin{pmatrix} {\mathrm{e}}^{{\mathrm{j}}\pi \tau_{\text{pmd}}  f } & 0 \\ 0 & {\mathrm{e}}^{-{\mathrm{j}}\pi \tau_{\text{pmd}}  f } \end{pmatrix}  {\mathrm{e}}^{-\j2\pi^2 \beta_{\text{cd}} L_{\text{cd}}  f^2 }  , \qquad \text{with} \quad \bm{R} = \begin{pmatrix} \cos\mlr{\gamma_{\text{hv}}} & \sin\mlr{\gamma_{\text{hv}}} \\ -\sin\mlr{\gamma_{\text{hv}}} & \cos\mlr{\gamma_{\text{hv}}} \end{pmatrix} \ ,
\end{align*}
where $\bm{R}$ represents a static rotation of the reference polarization to the fiber's \ac{PSP}, called HV-phase-shift $\gamma_{\text{hv}}$. In contrast to \cite{lauinger2022blind}, we do not consider the initial HV-phase-shift as compensated yet.
Additionally, we include first-order \ac{PMD} caused by the differential group delay $\tau_{\mathrm{pmd}}$ between the \acp{PSP},
and residual \ac{CD}, which is defined by the fiber's group velocity dispersion (GVD) parameter $\beta_{\text{cd}}=-\SI{26}{\pico\second^2\per\kilo\meter}$ (equals $D_{\text{cd}} = \SI{20}{\pico\s\per\nano\m\per\kilo\m}$ at $\lambda=\SI{1550}{\nano\m}$) times the uncompensated fiber length $L_{\mathrm{cd}}$.
Complex \Ac{AWGN} is added on both polarizations.

We evaluate $N_\text{ind}=100$~frames ($N_\text{ind}=20$ for the VAEflex/CMAflex schemes)
and apply the same averaging scheme as in \cite{lauinger2022blind}, which results in an averaging over $80\,000\ldots100\,000$ symbols per frame index and run.
As performance metric, we use both the \ac{BMI} (also often denotes as GMI) and  the percentage of \emph{failed} runs, which we define as runs where $\text{BMI}<5$ ($\text{BMI}<4.8$ for \ac{PCS} with $\mathcal{H}=5.73$, respectively) after a sufficient number of samples.
Similarly to~\cite{lauinger2022blind}, we apply a scheduler which halves the learning rate after every 20 (5 for the VAEflex/CMAflex schemes) frames for all algorithms.

\begin{figure} [!t]
		\centering
		\includegraphics{tikz_ext/pol_full.pdf}
		\vspace*{-3ex}
		\caption{Results for different initial HV-shift $\gamma_{\mathrm{hv}}$ (left), and the convergence behavior at $\gamma_{\mathrm{hv}}=\num{0.25}\pi$ (right), where the solid (-----) curve corresponds to the convergence of a typical run, the dotted ($\mathbf{\cdot\!\cdot\!\cdot\!\cdot\!\cdot\!}$) to the convergence of a bad run, and the dash-dotted ($\mathbf{\cdot-\cdot}$) to the convergence of a good run. }
		\label{Results:pol}
		\vspace*{-4ex}
\end{figure}
\vspace*{-1ex}
\section{Numerical Results}\vspace*{-1.5ex}
All numerical results are simulated with $M_{\mathrm{eq}}=15$ and with a symbol rate of $R_{\text{S}}=\SI{100}{\giga Bd}$, an oversampling factor of $N_\text{os}=2$~\ac{sps}, and an SNR of $\SI{24}{dB}$ for uniform and $\SI{22}{dB}$ for \ac{PCS}-$64$-QAM with $\mathcal{H}=5.73$.

First, we analyze the convergence behavior at a challenging working point with $\gamma_{\mathrm{hv}}=\num{0.2}\pi$, $\tau_{\mathrm{pmd}} = \frac{T_\text{S}}2 = \frac1{2 R_\text{S}}$, and for various lengths $L_{\mathrm{cd}}$ of uncompensated fiber as depicted in Fig.~\ref{Results:Lcd_full}. To take the potentially long channel impulse response into account, we use $M_{\mathrm{est}}=25$~taps.
While the \ac{VAE} equalizers achieve the highest \ac{BMI}, the CMAbatch is the most stable algorithm for uniform QAM, where it converges reliably and within 5 frames as depicted in the right plot.
The CMAflex/VAEflex converge as fast as the CMAbatch, but we observe a higher percentage of failed runs (left plots) for higher \ac{CD}. However, the VAEflex achieves a significantly higher BMI than the CMAbatch and CMAflex. The computationally less expensive VAEbatch algorithm converges similarly reliably and equalizes well, as indicated by the constellation diagram in Fig.~\ref{Results:Lcd_full} as well as the high BMI.
The bottom row depicts the results for moderately shaped \ac{PCS}-64-QAM with an entropy of $\mathcal{H}=5.73$, which is commonly-known as challenging task for the popular \ac{CMA} based blind equalizers. Only the VAE based equalizers, which consider the \textit{a~priori} density, converge reliably until the \ac{CD} becomes too strong, while the CMAflex fails completely. The VAEflex is still converging within 10 frames but the VAEbatch converges more reliably.

At $\gamma_{\mathrm{hv}}=\frac{\pi}4$, the \ac{CMA}'s loss function has a local minimum at an invalid constellation, so the left plot of Fig.~\ref{Results:pol} depicts a sweep of $\gamma_{\mathrm{hv}}$ around that point. We set a rather strong \ac{PMD} with $\tau_{\mathrm{pmd}}=T_\text{S}$, a moderate $L_{\mathrm{cd}}=\SI{1}{\kilo\meter}$, and use $M_{\mathrm{est}}=15$. While the BMI after convergence stays similar (dashed), the percentage of failed runs (solid) increases significantly at $\gamma_{\mathrm{hv}}\approx\frac{\pi}4$  for all \ac{CMA} based equalizers, which eventually fail in over 50\% of the runs. However, the VAE based equalizers do not struggle at all for this working point.
The right plot shows the convergence behavior at the critical point of $\gamma_{\mathrm{hv}}=\frac{\pi}4$. We chose three runs per algorithm representing a good (dash-dotted) and a bad (dotted) performance as well as typical curve (solid) of the converged runs. While the \ac{CMA} based equalizers quickly reach a moderate BMI and converge slowly but gradually, the VAE based algorithms start at low BMI but converge rapidly. Hence, they reach their maximum BMI faster as all CMA based equalizers.
The results show that all algorithms are significantly affected by the learning rate scheduler, which might be subject of further optimization.

\vspace*{-1ex}
\section{Conclusion}\vspace*{-1.5ex}
We reveal critical working points in which the widely-used \ac{CMA} and its variants fail to converge. The novel \ac{VAE} based equalizers converge reliably at these critical points with both standard and \ac{PCS} formats.
The novel VAEflex is a fast-converging option and, for the CMA, we found that batch-wise updating is advantageous.

\vspace*{2ex}

\footnotesize
\noindent{\bf{Acknowledgements}}:
This work was carried out in the framework of
the CELTIC-NEXT project AI-NET-ANTILLAS (C2019/3-3) and was funded
by the German Federal Ministry of Education and Research (BMBF) under
grant agreement 16KIS1316.

\vspace*{-1ex}


\section{Conclusion}
\label{conclusion}

In this work,  a new data-driven approach to generating a Koopman  linear model based on the direct encoding of Koopman Operator (DDE) is presented as an alternative to dynamic mode decomposition (DMD) and other related methods using least squares estimate (LSE). The major contributions include: 1) The analytical formula of Direct Encoding is converted to a numerical formula for computing the inner product integrals from given data; 2) An efficient algorithm is developed and its convergence conditions to the true results are analyzed; 3) Numerical experiments demonstrate a) greater accuracy compared to EDMD, b) lower sensitivity to data distribution, and c) rapid convergence of inner product computation. Furthermore, the DDE method is incorporated to Deep Koopman, i.e. neural network based methods for construction of the Koopman Operator, for improving prediction accuracy.

\section{Experiments}
\label{experiments}

In this section, the DDE algorithm is implemented for the sake of evaluating its validity and comparing its modeling accuracy to EDMD.
Consider a 2nd order nonlinear system consisting of a pendulum with a nonlinear damper. See Fig. \ref{fig:penddiagram}. The pendulum also bounces against  walls with nonlinear compliance. The state variables for this system are $x = [ \theta, \dot{\theta}]^T$, and the equation of motion can be written as:
\begin{equation}
    \ddot{\theta} = - sin(\theta) + F_k + F_c
\end{equation}
where $F_k$ and $F_c$ are wall reaction moment and damping moment, respectively,
\begin{align}
    F_k &= \begin{cases}-\text{sign}(\theta) \: k(|\theta| - \frac{\pi}{4})^2 &\text{if }|\theta| \geq \frac{\pi}{4}\\
    0 &\text{if }|\theta| < \frac{\pi}{4}
    \end{cases}\\
    F_c &= -\text{sign}(\dot{\theta}) \: c \:\dot{\theta}^2
\end{align}
where $k = 200$ and $c = 1$.
We present a two part numerical experiment for this system. The first experiment regards variations in dataset size and distribution, and the second experiment varies the usage of observable functions.

\begin{figure}
    \centering
    \includegraphics[width = 0.9\linewidth]{Figures/Diagrams/pendulum.png}
    \caption{Diagram of pendulum with with walls. (a) depicts the range of the pendulum where the walls are equally angularly displaced from the vertical. (b) depicts the forces exerted on the pendulum due to contact with walls. (c) depicts the damping force exerted onto the pendulum. }
    \label{fig:penddiagram}
\end{figure}

\subsection{Dataset Variations}

The datasets tested are of three types:
\begin{enumerate}
\item \textbf{Uniform}: These datasets are composed of a rectangular dynamic range which is sampled uniformly, like an evenly divided grid. The range varies from $\theta = [-
0.8,0.8]$ and $\dot{\theta} = [-2,2]$, where the mass can hit the walls and the damping can vary from 0 to a significant value. See Fig. \ref{fig:penddiagram}-(b), (c).
\item \textbf{Gaussian}: Data points are sampled with a finite-support Gaussian distribution. The data are distributed non-uniformly with their highest density at the peak of the Gaussian placed at diverse locations. In addition, each dataset contains 100 data points uniformly distributed along the boarder of the dynamic range to guarantee the same dynamic range as the uniform datasets.  Samples outside the dynamic range are excluded.
\item \textbf{Trajectories}: These datasets are composed of trajectories, beginning from 100 initial conditions that are simulated forward the same number of time steps. The dynamic range of this dataset differs from the two other dataset types.
\end{enumerate}

The models constructed for DDE and EDMD use the same observable functions. The observable functions chosen are two dimensional radial basis functions (RBFs), uniformly distributed between the maximum and minimum values of each state variable in their respective dataset, and the state variables. The total order of the system is 27th order with 25 RBFs and 2 state variables.

A trajectory dataset graph is generated using Delaunay Triangles in DDE, shown in Fig. \ref{fig:pend_graph}.
\begin{figure}
    \centering
    \includegraphics[width= 0.9\linewidth]{Figures/pendulum/graph_shaded.png}
    \caption{Graph of connections for a trajectory dataset composed of 10000 points formed when utilizing the data-driven direct encoding method. The lightly red shaded region denotes the dynamic range. In the zoomed-in image, the difference in volumes associated to different data points can be observed. Noting the difference in size of the purple triangle versus the green triangle.}
    \label{fig:pend_graph}
\end{figure}
\begin{table}
\caption{Sum of Squared Errors over dynamic range with varying dataset sizes.}
\label{tab:pred_dynrange} \begin{center} \begin{tabular}{c c c}
\hline \hline

 \multirow{ 2}{*}{Dataset Size} & Total SSE & SSE Variance\\
 & EDMD / DDE  & EDMD / DDE\\
\hline
\hline
\multicolumn{3}{c }{\textit{Uniform Datasets}} \\
\hline
\hline
900  & 19.470 / \textbf{17.167}  &  \textbf{0.0094} / 0.0097\\
2500 & 17.995 / \textbf{16.471} &  \textbf{0.0095} / 0.0103\\
10000 & 17.010 / \textbf{16.184} & \textbf{0.0099} / 0.0105\\
22500 & 16.698 / \textbf{16.133}  &  \textbf{0.0101} / 0.0105\\
\hline
\hline
\multicolumn{3}{c}{\textit{Trajectory Datasets}} \\
\hline
\hline
 1000  & 56.532 / \textbf{33.392} & 0.0349 / \textbf{0.0193}\\
2500 & 33.330 / \textbf{25.064} &  0.0200 / \textbf{0.0130}\\
5000 & 31.690 / \textbf{25.099} &  0.0195 / \textbf{0.0129}\\
10000 & 30.184 / \textbf{25.101} &  0.0186 / \textbf{0.0129}\\
25000 & 29.380 / \textbf{25.106} &  0.0181 / \textbf{0.0129}\\
\hline
\hline
\end{tabular} \end{center} \end{table}

The accuracy of the models is tested through calculating sum of squared errors (SSE) for one-step ahead predictions over the dynamic range of the datasets. These error values are calculated for a uniform grid of points, similar to that used in the uniform datasets. A visualization of the SSE is plotted in Fig. \ref{fig:sse_space}. The results of these calculations are shown in Table \ref{tab:pred_dynrange} and \ref{tab:pred_gaussian}. In the computation, the dynamic range is discretized, and the SSE value of each point is summed. For the Gaussian datasets, the test is run for eight iterations of each dataset to account for randomness and the average result is noted.
\begin{figure}
    \centering
    \includegraphics[width = 1\linewidth]{Figures/pendulum/heatmapcomp.png}
    \caption{Sum of Squared Error plots for various datasets. (a) and (b) are EDMD models for Gaussian datasets; (a) uses a dataset that is centered at $[0,0]$ and (b) uses a dataset that is centered at $[0,2]$. (c) DDE and (d) EDMD models using the trajectory dataset composed of 2500 data points and using 25 RBFs. Circled in red are the regions with greatest variation in SSE between models. Only the dynamic range is shown in all plots.}
    \label{fig:sse_space}
\end{figure}
\begin{table}
\caption{Sum of Squared Errors over dynamic range for Gaussian distributions with different centers. Centers of the distribution are noted above each column.}
\label{tab:pred_gaussian} \begin{center} \begin{tabular}{c c c c}
\hline \hline
 \multirow{ 3}{*}{Dataset Size} & \multicolumn{3}{c }{Total SSE} \\
 & Center: $[0,0]$ & Center: $[0.8, 0]$ & Center: $[0,2]$ \\
 & EDMD / DDE  & EDMD / DDE & EDMD / DDE\\
\hline \hline
\multicolumn{4}{c}{\textit{Gaussian Datasets}} \\
\hline
\hline
1000  & 25.565 / \textbf{24.377} &  25.161 / \textbf{22.98} & 24.338 / \textbf{22.445}\\
2500 & 24.991 / \textbf{23.739}&  25.421 / \textbf{22.747} & 24.221 / \textbf{21.590}\\
 5000 & 24.662 / \textbf{23.518} &  26.805 / \textbf{22.415} & 23.914 / \textbf{21.195}\\
10000 & 24.395 / \textbf{23.085} &  28.424 / \textbf{21.825} & 25.770 / \textbf{21.052}\\
25000 & 25.219 / \textbf{22.167} &  30.788 / \textbf{21.687} & 26.941 / \textbf{20.704}\\
\hline \hline
\end{tabular} \end{center} \end{table}

\begin{table} [!ht]
\caption{Sum of Squared Errors over dynamic range with varying order of lifted linear models.}
\label{tab:pred_obs} \begin{center} \begin{tabular}{c c c}
\hline \hline
 \multirow{ 2}{*}{\# Observables} & \multicolumn{2}{c}{Total SSE}\\
  & EDMD & DDE \\
\hline
\hline
\multicolumn{3}{c}{\textit{Trajectory Dataset, 5000 points}} \\
\hline
27  & 31.690 & \textbf{25.099} \\
51  & 36.657 & \textbf{21.637}\\
83 & 28.437 & \textbf{13.613}  \\
\hline
\end{tabular} \end{center} \end{table}

\subsection{Observable Function Variation}

The second experiment varies the number of observable functions selected, thus increasing the order. In this experiment, the number of RBFs is varied through uniformly increasing the density of the centers of the function over the range of the dataset. The results are noted in Table \ref{tab:pred_obs}. In the table, the number detailing number of observables is the number of functions including the state variables.

\subsection{Discussion}
From these results we can make the following observations.

\begin{itemize}
    \item All the numerical experiments show that DDE outperforms EDMD in total SSE.
    \item For uniform datasets, both models are nearly equivalent, though DDE has slightly lower SSE in all cases, shown in Table \ref{tab:pred_dynrange}. This result is expected as all data points are weighted equally in a uniform distribution.
    \item EDMD models exhibit significantly different distributions of prediction error, depending on dataset distribution. In Fig. \ref{fig:sse_space}-(a), the EDMD model learned from a dataset with high density near the origin produces a prediction error distribution that is low in accuracy in the top-right and the bottom-left corners of the dynamic range. In contrast, Fig.\ref{fig:sse_space}-(b)  illustrates that when using EDMD to learn from data with high density at $ \theta = 0, \dot{\theta} = 2$, the model results in low accuracy in the lower half of the dynamic range. DDE does not exhibit this high distribution dependency and achieves lower total SSE, as shown in Table \ref{tab:pred_gaussian}.
    \item In the trajectory datasets in Table \ref{tab:pred_dynrange}, DDE is not only lower in total SSE than EDMD but is also smaller in variance. Fig. \ref{fig:sse_space}-(c) shows that DDE has a uniformly low error distribution across the dynamic range, while EDMD in Fig.\ref{fig:sse_space}-(d) has two regions, as circled in the figure, with significantly larger error.
    These regions coincide with sparsity in the dataset, providing evidence of EDMD's bias towards regions of high data density and explaining the difference in performance between models for the non-uniform datasets.
    \item In the trajectory datasets, the total SSE converges for DDE beginning with small dataset sizes. This result implies that the elements in the $R$ and $Q$ matrices of DDE, that is, the inner product integral computations, converge as the data size and the data density increase. This convergence is confirmed in Fig. \ref{fig:qplot}, where several elements of the $Q$ matrix are plotted against the data size.
    \item The second experiment, regarding variations in observable function numbers demonstrates the effect of DDE remains even for significant increases in observable functions, referring to Table \ref{tab:pred_obs}. In all cases tested, DDE significantly outperforms EDMD over the dynamic range, as expected for a non-uniform dataset.
\end{itemize}
\begin{figure}
    \centering
    \includegraphics[width = 0.85\linewidth]{Figures/pendulum/q_plot.png}
    \caption{Change in values of diagonal elements in the $Q$ matrix of DDE with respect to trajectory dataset size. The specific elements shown correspond to a state variable, $Q_{0,0}$, and three RBF functions of varying distances away from the equilibrium.}
    \label{fig:qplot}
\end{figure}

\section{Introduction}
\label{intro}

Dynamic Mode Decomposition (DMD) was presented as a method to produce linear models from data generated through nonlinear dynamical processes by using Singular Value Decomposition (SVD) \cite{schmid2010dynamic}. Later, this method was developed further to create Extended Dynamic Mode Decomposition (EDMD), which introduced the concept of using observable functions, nonlinear functions of state variables, as a method of augmenting the state space \cite{williams2015data}. EDMD referenced the Koopman Operator as justification and a theoretical underpinning for lifting the state space. In his seminal work, Bernard Koopman showed the existence of this operator that transforms nonlinear systems into linear systems \cite{koopman1931hamiltonian}. Another extension of DMD has shown the viability of using DMD for control on non-autonomous systems \cite{proctor2016dynamic}. This enabled complex nonlinear Model Predictive Control (MPC) to be converted to linear MPC \cite{Korda_2018}, leading to numerous studies  utilizing the methodology to real systems\cite{ng2021model, cibulkaVehicle2020, calderon2021, bruder2021}.

To improve the accuracy of the models based on the Koopman Operator, two avenues of research have formed. The first avenue regards the selection of the observable functions used for constructing a lifted state space, as these functions are a key ingredient in creating an accurate linear model. Various methods have been developed, including deep neural networks for learning effective observable functions \cite{Lusch_2018, yeung2017, Selby_2021} and optimization \cite{korda2020optimal}. However, an efficient selection of observables does not solve all the issues that arise when attempting to construct an accurate linear model. For example, it is known that unstable modes are involved in Koopman-based DMD models and their extensions although the underlying nonlinear systems are known to be stable. The second avenue of research involves the formulation of the linear transition matrix. Extensive studies have been done to create stable linear models to remedy the situations where an outright use of DMD would lead to the creation of an unstable linear model \cite{fanstabkoop2021, mamakoukas2020learning, bevanda2022diffeomorphically, han2021desko}. Recently, an extension of DMD, called Robust Dynamic Mode Decomposition (RDMD), utilizes statistical measures to suppress the effect of outliers on modeling the linear Koopman matrix\cite{Hossein_Abolmasoumi_2022}.

A fundamental difficulty in constructing a proper linear model is data dependency.
Least Squares Estimation (LSE), involved in all DMD based methods, often produces a significant bias due to the distribution of the dataset. To eliminate this dependency on data distributions, the current work takes an alternative approach to LSE.

Recently, a new formulation of the Koopman Operator, termed Koopman Direct Encoding (DE), was produced \cite{asada2022}. This method directly encodes the nonlinear function of state transition using observables as basis functions to obtain a Koopman linear model. Inner products of observable functions in composition with the nonlinear state transition function are used to construct the state transition matrix without use of LSE. While DE theoretically guarantees the exact linear model, it requires access to the nonlinear state equations, which are often not available in practical applications. The current work aims to fill the gap between DE and data-driven approaches.

There are four significant contributions presented in this work. The first is the conversion of the DE formula of the Koopman Operator to a data-driven formula. The second is a computational algorithm and proof of its convergence to the true inner products that constitute the DE formula. The third is numerical experiments that provide evidence that the proposed method, unlike EDMD, does not exhibit biases to data distribution, but can produce consistently higher accuracy compared to EDMD.
Finally, the DDE algorithm is utilized in modeling a high order nonlinear system in combination with deep learning.

\section{Data-Driven Koopman Encoding}
\label{methods}

The prevailing method for construction of the Koopman Operator, EDMD, is based on LSE and SVD. This method, however, cannot provide a consistent estimate; the result is highly dependent on the distribution of data as the Koopman Operator is being approximated \cite{ljung1998system}. This dependency on distribution occurs because a core assumption of LSE is that the model structure is correct; when this assumption is violated, LSE is unable to create an unbiased estimator \cite{freedman2009statistical}. As the Koopman operator is truncated in practical use, this assumption does not hold.

Non-uniform data distributions inevitably occur in practical applications. For a nonlinear dynamical system with a stable equilibrium, for example, data collected from experiments and/or simulation of the system tend to be dense in the vicinity of the equilibrium, as all trajectories that begin within a region of attraction converge to the equilibrium.
Because LSE applies equal weighting to all data points, the model is heavily tuned to the behavior of the densely populated region.

The Direct Encoding method described previously enables us to obtain the exact linear state transition matrix $A$ through inner product computations. As the formulation is based on integration over the entire state space, there is no bias towards particular parts of the domain.

However, the original form of the Direct Encoding method utilizes the nonlinear state equation, i.e. the self-map $f(x)$, to compute the inner products. In practical applications, such a nonlinear function is not always available; only data are available. The objective of this section is to establish a computational algorithm to obtain the $A$ matrix by  numerically computing the inner products, $\langle g_i, g_j \rangle, \langle g_i \circ f, g_j \rangle$, from a given set of data.

The method presented consists of three operations.
\begin{itemize}
    \item The integral of the inner products is reduced in range from the entire state space to the dynamic range encapsulated by the data.
    \item The dynamic range is discretized with data points.
    \item The inner product integral is reduced to a weighted summation of the integrand evaluated at each data point multiplied by the volume $\Delta v$ associated to each point.
\end{itemize}
 Naturally, if data are densely populated in a small region, the discretized integral interval is small and thereby the volume also becomes small. Similarly, the volume tends to be larger where the data are sparse. In the summation, the integrand evaluated at individual data points are "weighted" by the size of the volume. This numerical inner product calculation prevents overemphasis of clustered data.

\subsection{Inner Product Computation}

We present the data-driven encoding method (DDE) as an alternative data-driven method to DMD for calculating a finite order approximation of the Koopman Operator.
The objective of this method is to compute the matrices $R$ and $Q$ in (13) and (14) from data. This entails the computation of inner products:
\begin{align}
    \langle g_i, g_j \rangle = \int_X G_{ij}(\xi)  d\xi \\
    \langle g_i \circ f, g_j \rangle = \int_X F_{ij}(\xi)  d\xi
\end{align}
where
\begin{equation}
    G_{ij} (x) = g_i(x) \bar{g}_j(x),  F_{ij} (x) = g_i[f(x)] \bar{g}_j(x)
\end{equation}
are assumed to be Riemann Integrable; the functions are bounded and continuous \cite{davis2007methods}.

There are two data sets used for the inner product computation. The first data set is
\begin{equation}
    D_N = \{ x_i \mid i = 1, \cdots, N; x_i \in X \}
\end{equation}
Note that all the data values are finite, $|x_i| < \infty$. As such, the integral interval of the inner products is finite in computing them from the data. To define the integral interval, we consider the dynamic range of the system, $X_D$, determined from the data set $D_N$. See Fig.\ref{fig:convexhull}. The dynamic range $X_D$ is defined to be the minimum domain in the space $X$ that includes all the data points in $D_N$, $X_D \supset D_N$, and that is convex. Namely, for any two states in $X_D$, $x, x' \in X_D$,
\begin{equation}
    \xi = \alpha x + (1-\alpha) x' \in X_D
\end{equation}
where $0 \leq \alpha \leq 1$.
\begin{figure} [!h]
    \centering
    \includegraphics[width = 0.8\linewidth]{Figures/Diagrams/convexhull.PNG}
    \caption{Illustration of the dynamic range of a dataset defined by the convex hull containing all points in the set, partitioned using a triangulation method. The data points are in black and the convex hull that encapsulates all data points is in grey.}
    \label{fig:convexhull}
\end{figure}
Each data point $x_i$ is mapped to $f(x_i)$, following the state transition law in eq.(\ref{eq:gen_nonlinear}). We assume that the transferred state, too, stays within the same dynamic range $X_D$. Collecting all the transferred states yields the second data set.
\begin{align}
    D_N^f = \{ f(x_i) \mid i = 1, \cdots, N; x_i \in D_N \} \\
    D_N^f \subset X_D
\end{align}
This implies that the state space of the nonlinear system under consideration is closed within the dynamic range $X_D$.

With this dynamic range, we redefine our objective to compute the inner products over $X_D$.
\begin{align}
    R_{ij} = \int_{X_D} G_{ij}(x) dx \label{eq:G_ij integral} \\
    Q_{ij} = \int_{X_D} F_{ij}(x) dx \label{eq:F_ij integral}
\end{align}
The integrals can be computed by partitioning the domain $X_D$ into many segments $X_1, \cdots X_P$, as shown in Fig. \ref{fig:convexhull}.
\begin{equation}
    X_D = \bigcup_{p=1}^P X_p
\end{equation}
We generate these segments by applying a meshing technique to the data set $D_N$, where the $n$-dimensional coordinates of individual data points are treated as nodes of a mesh. Delaunay Triangulation, for example, generates a triangular mesh structure with desirable properties \cite{cheng2013delaunay}. As illustrated in Fig. \ref{fig:convexhull}, each triangular element is convex and has no internal node. The volume of the dynamic range $V(X_D)$ is the sum of the volumes of all the elements.
\begin{equation}
    V(X_D) = \sum_{p=1}^P \Delta v_p
\end{equation}
Accordingly, the integral $R_{ij}$ in eq.(\ref{eq:G_ij integral}) can be segmented to
\begin{equation}
    R_{ij} = \sum_{p=1}^P \int_{X_p} G_{ij} (x) dx
\end{equation}
Suppose that the $p$-th element has $K_p$ nodes, as shown in Fig.\ref{fig:partition}. Renumbering these nodes 1 through $K_p$,
\begin{equation}
    \{ x[k_p] \Big| x[k_p] \in X_p; k_p = 1, \cdots, K_p \},  p = 1, \cdots, P
\end{equation}
The integrand $G_{ij}$ within the $p$-th element can be approximated to the mean of the $K_p$ nodes involved in the $p$-th element.
\begin{equation}
    G_{ij} (x;p) \approx \bar{G}_{ij,p}= \frac{1}{K_p} \sum_{k_p = 1}^{K_p} G_{ij}(x[k_p]), x[k_p] \in X_p
    \label{integrand_approx}
\end{equation}
If Delaunay Triangulation is used, $K_p = n+1$.
See Fig. \ref{fig:partition}. Substituting this into (\ref{eq:G_ij integral}) yields the approximate value of $R_{ij}$.
\begin{equation}
    \hat{R}_{ij} = \sum_{p=1}^P \frac{1}{K_p} \sum_{k_p = 1}^{K_p} G_{ij}(x[k_p])  \Delta v_p
    \label{eq:r_calc}
\end{equation}
where
\begin{equation}
     \Delta v_p = \int_{X_p} 1 \cdot dx
\end{equation}

Similarly, each component of the matrix $Q$ can be computed by using the same meshing.
\begin{equation}
    \hat{Q}_{ij} = \sum_{p=1}^P \frac{1}{K_p} \sum_{k_p = 1}^{K_p} F_{ij}(x[k_p]) \Delta v_p
    \label{eq:q_calc}
\end{equation}
Note that $F_{ij}$ is evaluated by using the data points in both $D_N^f$ and $D_N$,
\begin{equation}
    F_{ij}(x[k_p]) = g_i[f(x[k_p])] \bar{g}_j(x[k_p])
\end{equation}
where $f(x[k_p]) \in D_N^f$.
\begin{figure}
    \centering
    \includegraphics[width = 0.8\linewidth]{Figures/Diagrams/partition.PNG}
    \caption{Visualization of the integrand calculation process. The volume of the partition is encapsulated by the data points is denoted as $\Delta v_p$. With this grouping, the value of $G_{ij}$ is calculated for each point and the average among this group is computed, $\bar{G}_{ij,p}$.  In turn, this value, weighted by the volume of this partition, is summed across other partitions (not shown) to approximate the value of the element $R_{ij}$.}
    \label{fig:partition}
\end{figure}

\subsection {Convergence}

Consider the center of each partition, $\bar{x}_p = \int_{X_p} x dx / \Delta v_p$, and the distance between $\bar{x}_p$ and each point, $x[k_p]$:
\begin{equation}
    \Delta x[k_p] = \bar{x}_p - x[k_p]
\end{equation}
See Fig.\ref{fig:partition}. The maximum distance from the center of the partition to each point that makes up the partition is
\begin{equation}
    |\Delta x_p| = \max\{|\Delta x[1]|,\, ... \, |\Delta x[k_p-1]|, \, |\Delta x[k_p]| \}
\end{equation}
Consider a sequence of refining the approximate inner product integral $\hat{R}_{ij}$ by increasing data points $N$.  We can show that, as the number of partition $P$ tends infinity and the maximum subintervals $|\Delta x_p|$ approach zero, the approximate inner product integral $\hat{R}_{ij}$ converges to its true integral.

\begin{equation}
    R_{ij} =  \lim_{\substack{ P \to \infty \\ |\Delta x_p| \to 0}} \ \sum_{p=1}^P \frac{1}{K_p} \sum_{k_p = 1}^{K_p} G_{ij}(x[k_p])  \Delta v_p
\end{equation}

This formulation takes the form of weighted sums, specifically Riemann sums.
Given functions that are bounded and continuous over the subdomain of interest, sequences of this form are known to have a common limit and thus converge upon refinement to the Riemann integral value over that subdomain, according to Numerical Integration theory \cite[Section 1.5]{davis2007methods}.

\subsection{Algorithm}

In the prior section, integrals (\ref{eq:r_calc}) and (\ref{eq:q_calc}) are presented as summations over partitions. This computation can be streamlined by converting the summations over partitions to the one over nodes. Consider node 3 associated to data point $x_3$ in Fig.\ref{fig:convexhull}, for example. This node is an apex of the 5 surrounding triangles. This implies that integrand $G_{ij}(x)$ is calculated or recalled 5 times in computing (\ref{eq:r_calc}) and (\ref{eq:q_calc}). This repetition can be eliminated by computing volume $\Delta v_k$ associated to node $k$, rather than partition $p: \Delta v_p$. Namely, we compute
\begin{equation}
    \Delta v_k = \sum_{p=1}^{P} \frac{\Delta v_p}{K_p}I(k,p)
\end{equation}
where $I(k,p)$ is a membership function that takes value 1 when node $k$ is an apex of triangle $p$, that is, node $k$ is involved in partition $p$. Using this volume as a new weight we can rewrite (\ref{eq:r_calc}) and (\ref{eq:q_calc}) to be
\begin{align}
    \hat{R}_{ij} = \sum_{k=1}^K G_{ij}(x[k])  \Delta v_k  \label{eq:r_calc2}\\
    \hat{Q}_{ij} = \sum_{k=1}^K F_{ij}(x[k]) \Delta v_k  \label{eq:q_calc2}
\end{align}

Using this conversion, the computation can be streamlined and cleanly separated into three steps, as shown by pseudo-code in Algorithm \ref{algo:dde}.
The steps are: (1) \textit{Graph Creation:} data are connected to create partitions of the domain using a mesh generator: lines 3 to 8, (2) \textit{Weighting Calculation:} calculation of the weights for each data point: lines 10 to 17, and (3) \textit{Matrix Calculation:} the calculation of the $R$ and $Q$ matrices to find the  matrix $A$, lines 19 to 21.

\begin{algorithm} [!b]
    \caption{Algorithm for DDE in pseudocode}
    \label{algo:dde}
    \KwIn{}
    $D_N$, $D_N^f$\\
    \KwOut{}
    $P$: $p$, $K$: $k$, $R$, $Q$ $A$\\
    \textit{Graph Creation:} \\
    \For{$x_i$ in $D_N$ and corresponding $f(x_i)$ in $D_N^f$}{
    Create node, $k$\\
    Assign node attributes for current data point $k.x_t = x_t$ and $k.x_{t+1} = x_{t+1}$  \\
    Append node to node list $K$.
    }
    \textbf{then}{ \\
    \textit{Weighting Calculation:} \\
    Create list of triangles, $P$ using Delaunay Triangulation on node list $K$ using $k.{x_t}$}\\
    \For{$p$ in $P$}{
    Calculate volume, $V$, of $p$\\
    \For{$k$ corresponding to $K_p$}{
    Update volume of each node $k.\Delta v_k = k.\Delta v_k + \frac{V}{n+1}$}
    }
    \textbf{then}{\\
    \textit{Matrix Calculation:}\\
    Find $R$ and $Q$ via (\ref{eq:r_calc2}) and (\ref{eq:q_calc2})\\
    Find $A = Q R^{-1}$
    }
\end{algorithm}

\section{Application to Neural Net Koopman Modeling}
\label{nnkoop}
The use of deep neural networks for finding effective observable functions and constructing a Koopman linear model has been reported by several groups \cite{Lusch_2018, yeung2017, Selby_2021}. This method, sometimes referred to as Deep Koopman, is effective for approximating the Koopman operator to a low-order model, compared to the use of locally activated functions, such as RBFs, which scale poorly for high-order nonlinear systems. The proposed DDE method can be incorporated into Deep Koopman, further improving approximation accuracy.

Fig. \ref{fig:nn} shows the architecture of the neural network similar to the prior works \cite{Lusch_2018, yeung2017, Selby_2021}. The input layer receives training data of independent state variables. The successive hidden layers produce observable functions; these functions feed into the output layer consisting of linear activation functions. This linear output layer corresponds to the $A$ matrix that maps the observables of the current time to those of the next time step, i.e. the state transition in the lifted space. In the Deep Koopman approach, the output layer, that is, the $A$ matrix, is trained together with observable functions in the hidden layers. Using the observable functions learned from deep learning, the output layer weights are replaced with the $A$ matrix obtained from DDE, captioned in Fig. \ref{fig:nn}.

\begin{figure}
    \centering
    \includegraphics[width = 0.8\linewidth]{Figures/Diagrams/nn.PNG}
    \caption{Feedforward neural network model used to generate observable functions. The final layer of the neural network is a linear layer. In the standard Deep Koopman model this remains the same after the model is fully trained. However, with the DDE model, this final layer is recalculated using DDE by taking in the dataset used to train the model.}
    \label{fig:nn}
\end{figure}

The effectiveness of this Deep Koopman-DDE method is applied to a multi-cable manipulation system \cite{ng2021model}. A simplified single cable version is utilized consisting of six independent state variables. The network is constructed using PyTorch, trained with an Adam optimizer, to generate 40 observable functions. These model concatenates the state variables of the nonlinear system, resulting in a 46th order model. Table \ref{tab:nndef} shows parameters used for training the Deep Koopman model. The training dataset is composed of 3000 data points drawn from trajectories. Table \ref{table:SSE of Neural Net} compares the Deep Koopman model to the proposed model that uses DDE. Results are in terms of sum of squared error over a set of test trajectories. A significant improvement is achieved by incorporating DDE into the deep learning method.
\begin{table} [!h]
    \caption{Neural network parameters and characteristics.}
    \small
    \centering
    \begin{tabular}{c|c}
    \hline
        Parameters and Characteristics & Value \\
        \hline
        Number of Hidden layers & 3 \\
        Activation Functions, Both Hidden Layers & ReLU \\
         Width of 1st Hidden Layer & 16 \\
         Width of 2nd Hidden Layer & 16 \\
         Width of 3rd Hidden Layer & 40 \\
         Learning Rate, $\alpha$ & 0.01 \\
    \end{tabular}
    \label{tab:nndef}
\end{table}
\begin{table} [!h]
\caption{Average SSE prediction error for trajectories in set of test data for single winch system.}
\label{table:SSE of Neural Net} \begin{center} \begin{tabular}{c c c}
\hline \hline
Modeling Method & 1 Time Step & 20 Time Steps  \\
\hline
\hline
Deep Koopman only  & 0.2471 & 9.8610  \\
Deep Koopman + DDE  & \textbf{0.2350} & \textbf{4.1131}  \\
\hline
\end{tabular} \end{center} \end{table}

Practical concerns arose from the use of Delaunay Triangulation in this experiment. Specifically, the method was found to be limiting and inconsistent for systems that are eighth order or higher, due to computational cost and numerical instability. Alternative numerical integration approaches or alternative partitioning methods may be necessary when applying the method to higher order nonlinear systems.

\section{Koopman Operator and the Direct Encoding Method}
\label{preliminaries}

In this section, we give a brief overview of the Koopman Operator and dynamic mode decomposition, and introduce the direct encoding method for obtaining a Koopman operator directly from nonlinear dynamics.

\subsection{Least Squares Estimation of the Koopman operator}

Consider a discrete-time dynamical system, given by

\begin{equation}
    x_{t+1} = f(x_t)
    \label{eq:gen_nonlinear}
\end{equation}
where $x \in X \subset \mathbb{R}^n$ is the independent state variable vector representing the dynamic state of the system,
$f$ is a self-map, nonlinear function $f : X \rightarrow X$, and $t$ is the current time step. Also consider a real-valued observable function of the state variables $g : X \rightarrow \mathbb{R}$. The Koopman Operator $\mathcal{K}$ is an infinite-dimensional linear operator acting on the observable function $g$ :
\begin{equation}
    \mathcal{K} g  = g \circ f
\end{equation}
where $g \circ f$ is the composition of function $g$ with function $f$: $(g \circ f)(x) = g(f(x))$.

A common data-driven method for constructing the operator is Extended Dynamic Mode Decomposition (EDMD) \cite{williams2015data}, where observables that are experimentally obtained or simulated from the governing equation of the system are augmented by including real-valued observable functions of the independent state vector $x_t$. Collectively, a high-dimensional state vector is formed.
\begin{equation}
    z_t = \begin{bmatrix} g_1(x_t) \\ g_2(x_t) \\ \vdots \\ g_m(x_t) \end{bmatrix}
\end{equation}
where $m$ is the order of the lifted state corresponding to the number of observable functions. Underpinned by the Koopman Operator theory, EDMD assumes the existence of a linear state transition matrix $A$ relating $z_{t+1}$ to $z_t$, and determine $A$ by solving a least squares regression that minimizes the Sum of Squared Error (SSE).

\begin{equation}
   A  = \arg \min_A \sum_{t}|| z_{t+1} - A z_t || ^2
\end{equation}
Singular Value Decomposition (SVD) is used for the least squares optimization.
\subsection{Direct Encoding of the Koopman Operator}
\label{koop:de}

An alternative to the least squares estimate and EDMD is to obtain the exact $A$ matrix by directly encoding the self-map, nonlinear state transition function $f(x)$ with an independent and complete set of observable functions through inner product computations. This Direct Encoding method is introduced next, while the full proof can be found in \cite{asada2022}.

Let us first consider the case where $g_1, g_2, g_3, ...$ are orthonormal basis functions spanning a Hilbert space $\mathcal{H}$.
We assume that the self-map nonlinear function $f(x)$ is continuous and that the composition of $g_j$ with $f$ is also involved in the Hilbert space.
\begin{equation}
    g_j \circ f \in \mathcal{H} \quad\forall j
\end{equation}
This implies that the function $g_j \circ f$ can be expanded in $[g_1, g_2, g_3, ...]$.
\begin{equation}
    g_j \circ f = \sum_{k=1}^\infty \langle g_j \circ f, g_k \rangle g_k
    \label{eq:orthof_exp}
\end{equation}

Concatenating $g_1, g_2, g_3, ...$ and $g_1 \circ f, g_2 \circ f, g_3 \circ f ...$ in infinite dimensional column vectors, respectively,
\begin{equation}
   \bar{z}_t = \begin{bmatrix} g_1(x_t) \\ g_2(x_t) \\ \vdots  \end{bmatrix},
   \bar{z}_{t+1} = \begin{bmatrix} g_1[f(x_t)] \\ g_2[f(x_t)] \\ \vdots  \end{bmatrix}
\end{equation}
eq. (\ref{eq:orthof_exp}) can be written in matrix and vector form.
\begin{equation}
    \bar{z}_{t+1} = \bar{A} \bar{z}_t
    \label{eq:X_t+1 = A X_t}
\end{equation}
where $\bar{A}$ is an infinite dimensional matrix consisting of the inner products involved in eq. (\ref{eq:orthof_exp}),

\begin{equation}
    \bar{A} = \begin{bmatrix} \langle g_1\circ f, g_1 \rangle & \langle g_1 \circ f, g_2 \rangle & \dots \\ \langle g_2 \circ f, g_1 \rangle & \langle g_2 \circ f, g_2 \rangle & \dots \\ \vdots & \vdots & \ddots \end{bmatrix}
\end{equation}
Eq. (\ref{eq:X_t+1 = A X_t}) manifests that the state lifted to the infinite dimensional space $\bar{z_t}$ makes linear state transition with matrix $\bar{A}$.

The observables $g_1, g_2, g_3, ...$ were assumed to be orthonormal basis functions in the above derivation. This assumption can be relaxed to an independent and complete set of basis functions spanning the Hilbert space. Hereafter, let $[g_1, g_2, g_3, ...]$ be an independent and complete set of basis functions spanning the Hilbert space.

It can be shown that the time evolution of lifted state $z_t$ is given by a constant matrix $A_f$ for the independent and complete set of basis functions $[g_1, g_2, g_3, ...]$.
\begin{equation}
    z_{t+1} = A_f z_t
    \label{eq:A_f 1}
\end{equation}

    \label{A_f with C}

The matrix $A_f$ can be computed directly from the self-map, state transition function $f(x)$ and an independent and complete set of observables $[g_1, g_2, g_3, ...]$ through inner product computations. Post-multiplying the transpose of $z_t$ to both sides of eq. (\ref{eq:A_f 1}) and integrating them over $X$ yield:
\begin{equation}
    \int_X z(f(x)) z^T(x)dx = A_f \int_X z(x) z^T(x) dx
\end{equation}
which can be written as
\begin{equation}
    Q = A_f R
\end{equation}
where
\begin{align}
    Q &= \begin{bmatrix} \langle g_1\circ f, g_1 \rangle & \langle g_1 \circ f, g_2 \rangle & \dots \\ \langle g_2 \circ f, g_1 \rangle & \langle g_2 \circ f, g_2 \rangle & \dots \\ \vdots & \vdots & \ddots \end{bmatrix}\\
    R &= \begin{bmatrix} \langle g_1, g_1 \rangle & \langle g_1, g_2 \rangle & \dots \\ \langle g_2, g_1 \rangle & \langle g_2, g_2 \rangle & \dots \\ \vdots & \vdots & \ddots \end{bmatrix}
\end{align}
Because the observables are independent, the matrix $R$ is non-singular. Therefore, the matrix $A_f$ is given by
\begin{equation}
    A_f = QR^{-1}
    \label{eq:A_f = Q*R_inv}
\end{equation}

This formula for obtaining the matrix $A_f$ directly from the governing nonlinear state equation with the function $f(x)$ and the independent observables through inner products, which are guaranteed to exist in Hilbert space $\mathcal{H}$, is the Direct Encoding method.