text:
\section{Introduction}

Codemixing is the phenomenon of intermixing linguistic units from two or more languages in a single utterance, and is especially widespread in multilingual societies across the world \cite{muysken2000bilingual}. With increasing internet access to such large populations of multilingual speakers, there is active ongoing research on processing codemixed texts on online socialmedia communities such as Twitter and Facebook \cite{singh-etal-2018-twitter,DBLP:AmeyaPrabhuJSV16}. Not only do these texts contain a diverse variety of language use spanning the formal and colloquial spectra, such texts also pose challenging problems such as out of vocabulary words, slangs, grammatical switching and structural inconsistencies.

Previously, various approaches \cite{DBLP:AmeyaPrabhuJSV16,DBLP:madanEmsemble,singh-etal-2018-twitter,singh2018twitter} have focused on task specific datasets and learning architectures for syntactic and semantic processing for codemixed texts. This has facilitated developement of various syntactic and semantic task specific datasets and neural architectures, but has been limited by the expensive efforts towards annotation. As a result, while these efforts have enabled processing of codemixed texts, they still suffer from data scarcity and poor representation learning, and the small individual dataset sizes usually limiting the model performance.

Curriculum Learning, as introduced by \cite{bengio2009curriculum} is \enquote{to start small, learn easier aspects of the task or easier subtasks, and then gradually increase the difficulty level}. They also draw parallels with human learning curriculum and education system, where different concepts are introduced in an order at different times, and has led to advancement in research towards animal training \cite{krueger2009flexible}. Previous experiments with tasks like language modelling \cite{bengio2009curriculum}, Dependency Parsing, and entailment \cite{DBLP:journals/corr/HashimotoXTS16} have shown faster convergence and performance gains by following a curriculum training regimen in the order of increasingly complicated syntactic and semantic tasks. \cite{DBLP:journals/corr/abs-1802-03796} also find theoretical and experimental evidence for curriculum learning by pretraining on another task leading to faster convergence.

With this purview, we propose a syntactico-semantic curriculum training strategy for Hi-En codemixed twitter sentiment analysis. We explore various pretraining strategies encompassing Language Identification, Part of Speech Tagging, and Language Modelling in different configurations. We investigate the role of different transfer learning strategies by changing learning rates and gradient freezing to prevent catastrophic forgetting and interference between source and target tasks. We also propose a new model for codemixed sentiment analysis based on character trigram sequences and pooling over time for representation learning. We investigate the convergence rate and model performance across various learning strategies, and find faster model convergence and performance gains on the test set.

\section{Related Work}
Research on semantic and syntactic processing of codemixed texts has increasingly gained attention, and various approaches have been proposed to this end. \cite{DBLP:AmeyaPrabhuJSV16} released a dataset comprising user comments on Facebook pages, and proposed a convolutions over character embeddings approach towards sentiment analysis for Hi-En Codemixed texts. \cite{DBLP:madanEmsemble} propose a character trigram approach coupled with an ensemble of an RNN and a Naive Bayes classifier towards sentiment analysis for codemixed data.

More generally for monolingual sentiment analysis, RNNs and other sequential deep learning models have shown to be successful. \cite{socher-etal-2012-semantic} obtained significant performance improvement by incorporating compositional vector representations over single vector representations. \cite{zheng2018leftcenterright} take a different approach by capturing the most important words on either side to perform targeted sentiment analysis. Their LSTM based model uses context2target attention to achieve better benchmark performance on three datasets.  \cite{singh2018twitter} developed a dataset for Hi-En codemixed Part of Speech tagging, and proposed a CRF based approach. \cite{singh2018twitter} developed a dataset for Hindi English Codemixed Language Identification and NER, and propose a CRF based approach with handcrafted features for Named Entity Recognition.

Bengio et. al.\shortcite{bengio2009curriculum} introduced curriculum learning approaches towards both vision and language related task, and show significant convergence and performance gains for language modelling task. \cite{DBLP:journals/corr/HashimotoXTS16} propose a hierarchical multitask neural architecture with the lower layers performing syntactic tasks, and the higher layers performing the more involved semantic tasks while using the lower layer predictions. \cite{DBLP:journals/corr/abs-1808-10485} also propose a syntactico semantic curriculum with chunking, semantic role labelling and coreference resolution, and show performance gains over strong baselines. Like \cite{DBLP:journals/corr/HashimotoXTS16}, they hypothesize the incorporation of simpler syntactic information into semantic tasks, and provide empirical evidence for the same.

\section{Datasets}

\cite{DBLP:AmeyaPrabhuJSV16} released a Hi-En codemixed dataset for Sentiment Analysis, comprising 3879 Facebook comments on public pages of Salman Khan and Narendra Modi. Comments are annotated as positive, negative and neutral based on their sentiment polarity, and are distributed across the 3 classes as 15\% negative, 50\% neutral and 35\% positive comments.

\cite{singh-etal-2018-twitter} released a twitter corpus for Part of Speech tagging for Hindi English codemixed tweets about 5 incidents, and annotated 1489 tweets with the POS tag for each token.

\cite{singh2018twitter} released a Hindi English codemixed twitter corpus for Language Identification and Named Entities, where each token is annotated for its Language as English, Hindi or rest. Apart from this, they also annotate 2763 named entities in text, classified as Person, Location or Organization.
\begin{table}[ht]
\centering
\begin{tabular}{@{}ll@{}}
\toprule
\textbf{Symbol}      & \textbf{Description} \\ \midrule
$h_{r_i}^{(j)}$ & jth Layer's right directional hidden state at time i    \\
$h_{l_i}^{(j)}$             &  jth Layer's left directional hidden state at time i \\
$H_{i}^{(j)}$ & Left and right hidden state concatenation \\
$H_{N}^{(j)} $ & Terminal H for layer j LSTM   \\
\bottomrule
\end{tabular}
\caption{Symbol Description}
\label{table:symbol}
\end{table}

\section{Approach}

In the following subsections, we introduce our proposed model architecture for processing the above tasks in a hierarchical manner, discuss the various curriculum strategies we experiment with, and finally discuss the transfer learning techniques we explore.

\begin{figure}[h]
\includegraphics[width=3.37in]{model}
\caption{Systematic Overview of Model Architecture}
\centering
\label{fig:model}
\end{figure}

\subsection{Model}

We case normalise the texts and mask user mentions and URLs with special characters. After tokenizing the texts, we append a token terminal \enquote{*} symbol to each token and further split each token into its constituent character trigrams. Thus, a token \enquote{girl} is split into \enquote{gir} + \enquote{l*\#}, where \enquote{\#} is the padding symbol for character trigrams. Due to the high class imbalance in Sentiment Analysis data, we perform a mixture of oversampling and undersampling between classes than simply prune the samples of the larger classes.

As shown in Figure \ref{fig:model}, the model comprises of an Embedding layer followed by two layers of bidirectional LSTMs. The embedding layer serves as a lookup table for our character trigram dense representations, and a sequence of these representations are passed on to the LSTM stack for each input sample.

Layer 1 of the LSTM stack takes the sequence of character trigram embeddings as input, and is used to predict the corresponding POS tag and the Language tag at each time step. The concatenation of the left and right directional hidden states at each time step i is passed to a standard softmax classifier to output the probability distribution over the POS tags at that timestep. Similarly, another softmax classifier outputs the language tags at each timestep.

\begin{flalign*}
&H_{i}^{(1)} =  concat([h_{r_i}^{(1)}; h_{l_i}^{(1)}])  &\\
&y_{i_{pos}} = Softmax(W_{pos}H_{i}^{(1)} + b_{pos}) &\\
&y_{i_{lang}} = Softmax(W_{lang}H_{i}^{(1)} + b_{lang}) &\\
\end{flalign*}

The LSTM Layer 2 takes the concatenated bidirectional hidden states of Layer 1 ($H_i^{(1)}$) as input to learn the sequence representation as an abstraction over the layer 1 representations. This architecture allows the semantic task to consider both the character trigram dense representation as well as "POS" and "Language" aware bidirectional representations to perform more complex semantic tasks like language modelling and sentiment analysis.

For language modelling, the concatenated terminal hidden states from the right and left directional LSTM layer 2 is passed to a standard softmax classifier, which outputs the probability distribution over the character trigram vocabulary for the next trigram in the input sequence. For sentiment analysis, we concatenate the max pooling over time, avg pooling over time, and the terminal hidden states of the Layer 2 BiLSTM to form the representation. The maxpooling and avgpooling over time representations circumvents the information loss in sequence terminal representations.  This representation is passed to a standard softmax classifier to predict the sentiment polarity over the 3 classes.

\begin{flalign*}
&H_{i}^{(2)} =  concat([h_{r_i}^{(2)}; h_{l_i}^{(2)}])  &\\
&y_{LM} = Softmax(W_{LM}H_{N}^{(2)} + b_{LM}) &\\
&H_{S}^{(2)} = concat([H_{N}^{(2)};maxpool(H_{i}^{(2)});avgpool(H_{i}^{(2)})]) &\\
&y_{sentiment} = Softmax(W_{sentiment}H_{S}^{(2)} + b_{sentiment}) &\\
\end{flalign*}

For each of the tasks described above, we train our model to optimize the cross entropy loss for the given prediction, formulated as :

\[L = -\sum_{y\in Y} y\log{p} + (1-y)log(1-p)\]

where y is the true label, and p is the predicted probability of that label by the model.

\subsection{Curriculum Training}

While our proposed model enables efficient transfer learning by progressive abstraction of representations for more complicated tasks, the highlight of the approach lies in the training regimen followed.

Curriculum learning can be seen as a sequence of training criteria \cite{bengio2009curriculum}, with increasing task or sample difficulty as the training progresses. It is also closely related with transfer learning by pretraining, especially in the case when the tasks form a logical hierarchy and contribute to the downstream tasks. With this purview, we propose a linguistic hierarchy of training tasks for codemixed languages, with further layers abstracting over the previous ones to achieve increasingly complicated tasks. Considering the codemixed nature of texts and linguistic hierarchy of information, we propose the tasks in the order of : Language Identification, Part of Speech Tagging, Language Modelling and further semantic tasks like sentiment analysis.

Since tokens in codemixed texts have distinct semantic spaces based on their source language, Language Identification can incorporate this disparity among the learnt trigram representations. Following this, the Part of Speech Tagging groups the words based on their logical semantic categories, and encodes simpler word category information in a sequence. Also, as in \cite{singh-etal-2018-twitter,sharma2016shallow}, Language Tag and Part of Speech Tag have previously been provided as manual handcrafted features for a range of downstream syntactic and semantic tasks.

In addition to the above tasks, Language Model pretraining has shown significant performance gains as reported by \cite{howard2018universal}. It captures various aspects of language such as long range dependencies \cite{linzen2016assessing}, word categories, and sentiment \cite{radford2017learning}.

Conforming with the linguistic hierarchical information, we first train our model to predict the language labels for each character trigram as per its token. This is followed by further training the model to predict the Part of Speech tag for each of its character trigram as per its token. Subsequently, the model is trained on Language Modelling task, in process training the LSTM Layer 2 to build over the LSTM Layer 1 inputs to learn meaningful sequence representation. Lastly, the model is trained to predict the sentiment of the input text based on the LSTM Layer 2 representation.

\begin{table}[t]
\centering
\begin{tabular}{@{}ll@{}}
\toprule
\textbf{Hyperparameter}      & \textbf{Value} \\ \midrule
Embedding Dimension & 64    \\
LSTM Cell Dimension & 64    \\
Dropout             & 0.2   \\
Learning Rate       & 0.04  \\
Batch Size          & 4     \\ \bottomrule
\end{tabular}
\caption{Model Hyperparameters}
\label{table:params}
\end{table}

\begin{table}[t]
\centering
\begin{tabular}{@{}ll@{}}
\toprule
\textbf{Character Encoding}      & \textbf{Accuracy} \\ \midrule
Character Unigram & 62.57    \\
Convolution over Character Unigrams & 64.23    \\
Character Trigram & 67.83    \\
Byte Pair Encoding(BPE)             & 64.61   \\ \bottomrule
\end{tabular}
\caption{Character Encoding Experiments}
\label{table:chenc}
\end{table}

{

\begin{table*}[ht]
\centering
\begin{tabular}{|l|c|c|c|c|}
\hline
\multicolumn{1}{|r|}{Model} & \multicolumn{1}{l|}{Accuracy} & \multicolumn{1}{l|}{Precision} & \multicolumn{1}{l|}{Recall} & \multicolumn{1}{l|}{F1-score} \\ \hline
SVM(Unigrams)               & 61.7                          & 0.579                          & 0.551                       & 0.565                         \\ \hline
SVM(Unigrams+Bigrams)       & 64.1                          & 0.609                          & 0.537                       & 0.566                         \\ \hline
MNB(Unigrams)               & 64.5                          & 0.748                          & 0.485                       & 0.588                         \\ \hline
MNB(Unigrams+Bigrams)       & 66.1                          & 0.698                          & 0.540                       & 0.609                         \\ \hline
SentiWordNet                & 51.5                          & -                              & -                           & 0.252                         \\ \hline
Char-trigram based LSTM \cite{DBLP:madanEmsemble}    & 65.2                          & 0.610                          & 0.563                       & 0.586                         \\ \hline
Vowel-Consonant based       & 62.8                          & 0.652                          & 0.522                       & 0.580                         \\ \hline
Sub-word\cite{DBLP:AmeyaPrabhuJSV16}                    & 69.2                          & 0.684                          & 0.623                       & 0.652                         \\ \hline
Our Approach                & 72.51\%          & 0.712           & 0.645        & 0.677          \\ \hline
\end{tabular}
\caption{Model Comparison}
\end{table*}

}

\begin{table}[t]
\centering
\begin{tabular}{@{}ll@{}}
\toprule
\textbf{Training}      & \textbf{Accuracy} \\ \midrule
From Scratch (No Curriculum) & 70.19    \\
+ POS + LangId Pretraining & 68.83    \\
+ LM Pretraining & 72.51    \\
LM only Pretraining & 72.16    \\
No Gradual Unfreezing             &  70.74  \\
No Discriminative Finetuning             & 71.68   \\
\bottomrule
\end{tabular}
\caption{Curriculum and Finetuning Experiments}
\label{table:curr}
\end{table}

\subsection{Transfer Learning}

As noted in earlier efforts \cite{howard2018universal} towards finetuning pretrained models for NLP tasks, aggressive finetuning can cause catastrophic forgetting, thus causing the model to simply fit over the target task and forget any capabilities gained during the pretraining stage. On the other hand, too cautious finetuning can cause slow convergence and overfitting. To this end, we experiment with different strategies which can be broadly categorized as:

\paragraph{Discriminative Finetuning :} As also noted by \cite{yosinski2014transferable}, different layers capture different types of information, and thus need to be optimised differently. In the context of our model, the embedding layer captures the individual character trigram information, the LSTM layer 1 is trained towards capturing the token level information such as Part of Speech and Language Tag, and the final LSTM Layer 2 is trained to capture the overall textual representation to perform Language Modelling and Sentiment Analysis. With this purview, similar to \cite{howard2018universal}, we propose optimizing different layers in our model to different extents, and keep lower step sizes for the deeper pretrained layers while finetuning on a downstream task. We thus split the parameters as $\{\theta_1, ..., \theta_l\}$ , where $\theta_i$ corresponds to the parameters of layer i, and optimize them with separate learning rates $\{\eta_1, ...., \eta_l\}$ . Also, when finetuning a pretrained layer for a downstream task, we keep $\eta_i < \eta_j; \forall i<j$.

 Thus, while finetuning the POS + Lang Id pretrained model for Language Modeling, we propose to keep the learning rates for Embedding Layer and LSTM Layer 1 lower than the LSTM Layer 2 weights. Similarly, when finetuning the Language Model for Sentiment Analysis, we keep the learning rates of the deeper layers lower than that of the shallower ones.

\paragraph{Gradual Unfreezing:} Similar to \cite{howard2018universal}, rather than updating all the layers together for finetuning, we explore gradual ordered unfreezing of layers. Thus, initially we freeze all the layers. Then starting from the last layer, we train the model for a certain number of epochs before unfreezing the layer below it. Thus for Sentiment Analysis finetuning, for the first epoch, only $\theta_{sentiment}$ receives the gradient updates, after which we unfreeze the $\theta_{lstm2}$, and subsequently unfreeze the lower layers in a similar manner.

\section{Experiments}

The input to the LSTM stack is the sequence of character trigram dense representations, which we keep as 64 dimensional vectors. We also explore other token representations such as sequence of unigrams, convolution over unigrams \cite{DBLP:AmeyaPrabhuJSV16}, and Byte Pair Encoding (BPE) \cite{sennrich2015neural}. BPE is an unsupervised approach towards subword decomposition, and has shown improvements in MT systems and summarization. We train our model from scratch for Sentiment Analysis using the above mentioned character encodings, and report the results in Table \ref{table:chenc}.

Our LSTM stack consists of two layers of bidirectional LSTMs, with 64 hidden state dimensions. We add a dropout layer with the dropout rate set to 0.2 between the LSTM layers to prevent overfitting. We experiment with average pooling and max pooling concatenation over hidden states for semantic prediction, similar to \cite{howard2018universal}, and observe increase in model accuracy by 2.2\% on sentiment analysis.

To evaluate our baseline for curriculum training experiments, we initially train the model from scratch on the single target task (Sentiment Analysis) for 25 epochs. We approach the evaluation of our curriculum by training the model sequentially for four subtasks - Language Identification, POS Tagging, Language Modelling and Sentiment Analysis. We evaluate the strategy of pretraining with only POS Tagging and Language Identification, and observe similar performance as no curriculum training. We hypothesize the potential reasons for this drop and find a significant divergence in character trigram occurance between the Source Tasks (POS + Lang Id) and Target Task(Sentiment Analysis). This experiement highlights the importance of inclusion of language model pretraining for better token level representation learning, and provides a better model prior for sequence representation (LSTM Layer 2 output). We experiment with only Language Modelling as pretraining task, and observe significant gains over no curriculum strategy.

We note the convergence of our model with and without curriculum training, and observe that the curriculum training regimen causes faster convergence, as has been observed in previous works \cite{bengio2009curriculum,howard2018universal}. This is expected as the model is pretrained on prior tasks already have a general purpose representation learning, and only needs to adapt to the idiosyncrasies of the target task, i.e. Sentiment Analysis in this case.

As discussed in Section 4.3, for our transfer learning optimization experiments, we segment the optimization of different parameters of our model with different learning rates, in order to limit catastrophic forgetting and interference among the tasks, as proposed by \cite{howard2018universal}. We segregate our model parameters in the following 4 groups:
\begin{itemize}
\item Emb Layer
\item LSTM Layer 1
\item LSTM Layer 2
\item Sentiment Linear Map
\end{itemize}
We set the learning rate of the previous layer $\eta_{l-1}$ = $\eta_l / 2.0$. For our gradual unfreezing experiments, we unfreeze the lower layer after training the model for 1 epoch with the lower layer unfrozen.

\section{Future Work}

In future work, we would like to explore word normalization and miltilingual embeddings in conjunction to learning representations from scratch. Another line of potential study could be investigation into why BPE is able to lead to performance gains as in monolingual domains, but fails in the codemixed multilingual tasks in our experiments. We would also like to explore convolution over character embeddings as a method to further circumvent the out of vocabulary problem with codemixed social media data.

We also plan to explore better representation learning for the semantic tasks. One particular direction we plan to explore  attention over the LSTM layer 2 as a weighted peek into the intermediate hidden states for the semantic classification task.

In future, we would like to experiment with more syntactic tasks like dependency label predictions, and also study more semantic tasks like aggression detection. Codemixed domains suffer from severe resource scarcity, and thus vocabulary divergence between various datasets proves as a roadblock to generalizable models, as observed in POS + Language Id pretraining experiments.

\bibliographystyle{named}


text:
\section{Supplementary Materials}

\subsection{S1. Bosonic Bogoliubov-de Gennes Hamiltonian}

The spin interaction Hamiltonian is given by
\begin{equation}\label{eqS1}
  H=\sum_{l}H_{l}+\sum_{i,l}[J_{\perp}+(-1)^{l+1}\delta J_{\perp}]\mathbf{S}_{i}^{l}\cdot\mathbf{S}_{i}^{l+1},
\end{equation}
and the intralayer spin Hamiltonian
\begin{eqnarray}
  H_{l}=&-&J_{\parallel}\sum_{\langle ij\rangle}\mathbf{S}_{i}^{l}\cdot\mathbf{S}_{j}^{l}
  +D\sum_{\langle\langle ij\rangle\rangle}\nu_{ij}\hat{z}\cdot(\mathbf{S}_{i}^{l}\times\mathbf{S}_{j}^{l}) \nonumber \\
  &+&\sum_{i}[K_{x}(S_{i,x}^{l})^2-K_{z}(S_{i,z}^{l})^2]. \label{eqS2}
\end{eqnarray}
Due to interlayer antiferromagntic (AFM) order, we apply the Holstein-Primakoff transformation $S_{i,+}^{1}=\sqrt{2S}a_{i}^{1}$, $S_{i,-}^{1}=\sqrt{2S}a_{i}^{1\dagger}$, $S_{i,z}^{1}=S-a_{i}^{1\dagger}a_{i}^{1}$ for first layer
and $S_{i,-}^{2}=\sqrt{2S}a_{i}^{2}$, $S_{i,+}^{2}=\sqrt{2S}a_{i}^{2\dagger}$, $S_{i,z}^{2}=-S+a_{i}^{2\dagger}a_{i}^{2}$ for second layer,
we can get the  Hamiltonian up to the quadratic term
\begin{equation}\label{eqS3}
  H=\sum_{i,l}[h_{0}a_{i}^{l\dagger}a_{i}^{l}+(K_{x}a_{i}^{l}a_{i}^{l}+h_{\perp}^{l}a_{i}^{l}a_{i}^{l+1}+H.c.)]-\sum_{\langle ij\rangle,l}J_{\parallel}S(a_{i}^{l\dagger}a_{j}^{l}+H.c.)+\sum_{\langle\langle ij\rangle\rangle,l}(-1)^{l}DS(\nu_{ij}a_{i}^{l\dagger}a_{j}^{l}+H.c.),
\end{equation}
where $h_{0}=(2K_{z}+K_{x}+3J_{\parallel}+2J_{\perp})S$ and $h_{\perp}^{l}=[J_{\perp}-(-1)^{l}\delta J_{\perp}]S$, $H.c.$ denote the Hermitian conjugate.
Considering the AFM interlayer and FM intralayer  coupling, we apply the Fourier transformation,
\begin{equation*}
  a_{i}^{l}=\frac{1}{N}\sum_{\mathbf{k}}e^{i\mathbf{k}_{\parallel}\cdot\mathbf{r}_{i}^{l}}e^{ik_{z}r_{i,z}^{l}}b_{\alpha,\mathbf{k}}^{l}, \quad a_{i}^{l\dagger}=\frac{1}{N}\sum_{\mathbf{k}}e^{-i\mathbf{k}_{\parallel}\cdot\mathbf{r}_{i}^{l}}e^{-ik_{z}r_{i,z}^{l}}b_{\alpha,\mathbf{k}}^{l\dagger};\quad  l\in odd
\end{equation*}
\begin{equation*}
  a_{i}^{l}=\frac{1}{N}\sum_{\mathbf{k}}e^{i\mathbf{k}_{\parallel}\cdot\mathbf{r}_{i}^{l}}e^{-ik_{z}r_{i,z}^{l}}b_{\alpha,\mathbf{k}}^{l}, \quad a_{i}^{l\dagger}=\frac{1}{N}\sum_{\mathbf{k}}e^{-i\mathbf{k}_{\parallel}\cdot\mathbf{r}_{i}^{l}}e^{ik_{z}r_{i,z}^{l}}b_{\alpha,\mathbf{k}}^{l\dagger}.\quad  l\in even
\end{equation*}
$\mathbf{k}_{\parallel}=(k_{x},k_{y})$. Then we can get the Hamiltonian in momentum space $H=\sum_{\mathbf{k}}\Psi_{\mathbf{k}}^{\dagger}H_{\mathbf{k}}\Psi_{\mathbf{k}}$ with the Bosonic BdG Hamiltonian
\begin{equation}\label{eqS4}
  H_{\mathbf{k}}=\left(\begin{array}{cccc}
  h_{1}(\mathbf{k}) & 0 & K_{x} & h_{\perp,k_{z}} \\
  0 & h_{2}(\mathbf{k}) & h_{\perp,k_{z}} & K_{x} \\
  K_{x} & h_{\perp,k_{z}}^{\dagger} & h_{1}^{T}(-\mathbf{k}) & 0 \\
  h_{\perp,k_{z}}^{\dagger} & K_{x} & 0 &  h_{2}^{T}(-\mathbf{k})
  \end{array}\right).
\end{equation}
where the basis is $\Psi_{\mathbf{k}}=(b_{A,\mathbf{k}}^{1},b_{B,\mathbf{k}}^{1},b_{A,\mathbf{k}}^{2},b_{B,\mathbf{k}}^{2},
b_{A,-\mathbf{k}}^{1\dagger},b_{B,-\mathbf{k}}^{1\dagger},b_{A,-\mathbf{k}}^{2\dagger},b_{B,-\mathbf{k}}^{2\dagger})^{T}$.
Note that $-\mathbf{k}=(-k_{x},-k_{y},k_{z})$ appears in the subscript of basis and in hole part of the Hamiltonian.
$h_{1}(\mathbf{k})=h_{0}\sigma_{0}+h_{x}\sigma_{x}+h_{y}\sigma_{y}+h_{z}\sigma_{z}$ and $h_{2}(\mathbf{k})=h_{0}\sigma_{0}+h_{x}\sigma_{x}+h_{y}\sigma_{y}-h_{z}\sigma_{z}$. $\sigma_{0}$ is the $2\times2$ identity matrice
and $\sigma_{i=x,y,z}$ are the Pauli matrices.
$h_{\perp,k_{z}}=J_{\perp}(1+e^{-ik_{z}})+\delta J_{\perp}(1-e^{-ik_{z}})$.

\subsection{S2. Calculation of the Z$_{2}$ number}
We first focus on the bilayer case. The Hamiltonian host a inversion symmetry: $\mathcal{P}=\tau_{z}s_{0}\sigma_{x}$. The eigenvalues
at four time-reversal-invariant momentum (TRIM), i.e. $\Gamma$ and three $M$ points (shown in Fig.~S1 (a)), are $1$, $1$, $-1$, $1$, respectively.
So that $(-1)^{Z_{z}}=-1$ indicate $Z_{1}=1$ for bilayer system.

Now we focus on the 3D case. Four Z2 invariants are needed to form
the combination [$\nu_{0}$,($\nu_{1}$,$\nu_{2}$,$\nu_{3}$)] to describe the topological phase in 3D and distinguish the "strong" and "weak".
Here it is hard to directly calculate the four topological invariant as there is no inversion symmetry. We enlarge the unit cell along the $y$ direction
and the Brillouin zone (BZ) is reduced to a rectangle in the $k_{x}$-$k_{y}$ plane.
We adopt the method introduced by Rui Yu, et al.~\cite{RYu} to calculate the evolution lines of Wannier centers
to get the topological invariant indirectly. The effective BZ adopted is shown in Fig.~S1 (b) as the shaded region.
At $k_{z}=0$ and $k_{z}=\pi$, the Wannier centers both winding one time, indicating $\nu_{0}=0$, $\nu_{3}=1$.
at $k_{x}=\pi$ and $k_{y}=\pi$, the Wannier centers both winding zero times, $\nu_{1}=0$, $\nu_{2}=0$.
Thus the bulk system is in a "weak" TI phase.

\begin{figure}
  \centering

  \includegraphics[width=0.7\textwidth]{figS1}\\
  \caption{(a). The 2D BZ and the reduced BZ. (b) The reduced 3D BZ. The shaded region is the region for computing the Z$_{2}$ number using the Wilson loop method. (c-f). The Wannier center evolution. at the four shaded region in (b).}\label{figs1}
\end{figure}

\begin{figure}[b]
  \centering

  \includegraphics[width=0.5\textwidth]{figS2}\\
  \caption{(a) Magnon bands of bilayer zigzag ribbons. The parameters are $J_{\perp}=0.3J_{\parallel}$, $D=0.05J_{\parallel}$, $K_{z}=0.2J_{\parallel}$,
  $K_{x}=J_{\parallel}$. In this and subsequent calculations, we adopted a compensated boundary.
  (b) Magnon eigenvalues in a parallelogram flake geometry. The inset shows the module square of wavefunction distribution of the corner modes.}\label{figs2}
\end{figure}

\subsection{S3. Corner Modes in bilayer}
The bilayer share a same Hamiltonian to Eq.~(3) in the main text just by setting $h_{0}=2K_{z}+K_{x}+3J_{\parallel}+J_{\perp}$
and $h_{\perp,k_{z}}=J_{\perp}$.
Then the bosonic Hamiltonian have an additional pseudo-inversion symmetry: $\mathcal{P}\tau_{z}H_{\mathbf{k}}\mathcal{P}=\tau_{z}H_{-\mathbf{k}}$
with $\mathcal{P}=\mathcal{P}^{-1}=\tau_{z}s_{0}\sigma_{x}$.
When $K_{x}=0$, the parity eigenvalues at $\Gamma$ and three M-points, i.e. the time-reversal invariant momentum (TRIM) points,
are $1$, $1$, $-1$, $1$, respectively, indicating $Z_{2}$ invariant $\nu=1$.
The bilayer system is thus topologically nontrivial and supports gapless edge states with two-components propagating oppositely at each edge.
The EPA breaks the pseudo-time-reversal symmetry and destroys the $Z_{2}$ topological phase and opens gap for zigzag-edged ribbon, as shown in Fig.~S2(a).
The same to the bulk case, the magnetic two-fold rotational symmetry $\mathcal{M}_{y}$ is preserved and rotation-graded topological invariant
along the line $k_{x}$ at $k_{y}=0$ is half-quantized, makes the system a second-order topological insulator.
The corner modes arise as a result, as shown in Fig.~S2(b) in a parallelogram flake.

\subsection{S4. The chiral Hinge modes at $K_{x}=0$}

In Fig.~S2 (a), we plot the magnon local density distribution of the chiral hinge modes. In (b), we plot the magnon bands along $k_{x}$-direction at $K_{x}=0$.
Without the EPA, there also exist the chiral hinge modes, coexisting with the gapless type-I TSMs, closing the gap.

\begin{figure}[t]
  \centering

  \includegraphics[width=0.5\textwidth]{figS3}\\
  \caption{(a). Magnon local density distribution at $k_{x}=\pi$ at upper surface hinges. (b). The magnon bands along $k_{x}$-direction at $K_{x}=0$.
  The other parameters are the same to Fig.~3 (b) in the main text.}\label{figS3}
\end{figure}

\begin{figure}[b]
  \centering

  \includegraphics[width=0.8\textwidth]{figS4}\\
  \caption{(a). The ribbon band of a bilayer honeycomb ferromagnetis with monoclinic stacking. The parameters are $J_{\perp}=0.2J_{\parallel}$,
  $D=0.05J_{\parallel}$, $K_{z}=0.2J_{\parallel}$, $K_{x}=J_{\parallel}$. (b). The corner modes in a parallelogram flake geometry. The inset
  shows the local density distribution of the down layer. The upper layers shares a similar distribution.
  (c). The magnon spectrum along $k_{z}$ direction with $\delta J_{\perp}=0.1J_{\parallel}$.}\label{figS4}
\end{figure}

\subsection{S6. Magnons in Monoclinic stacked honeycomb ferromagnets}

We first discuss the bilayer case to get an intuitive picture.
The monoclinic stacking of bilayer CrI$_{3}$ is a $\frac{1}{3}$-glide along the Zigzag direction compared with AA stacking.
Here we only consider the nearest interlayer AFM coupling.
The bosonic BdG Hamiltonian in $\mathbf{k}$-space is given by
\begin{equation}\label{eqS5}
  H_{\mathbf{k}}=\left(\begin{array}{cccc}
  h_{1}(\mathbf{k}) & 0 & K_{x} & h_{\perp,\mathbf{k}}^{0} \\
  0 & h_{2}(\mathbf{k}) & h_{\perp,-\mathbf{k}}^{0T} & K_{x} \\
  K_{x} & h_{\perp,-\mathbf{k}}^{0*} & h_{1}^{T}(-\mathbf{k}) & 0 \\
  h_{\perp,\mathbf{k}}^{0\dagger} & K_{x} & 0 &  h_{2}^{T}(-\mathbf{k})
  \end{array}\right),
\end{equation}
in the basis $\Psi_{\mathbf{k}}=(\psi_{\mathbf{k}},\psi_{-\mathbf{k}}^{\dagger})^{T}$ with
$\psi_{k}=(b_{A,\mathbf{k}}^{1},b_{B,\mathbf{k}}^{1},e^{i\frac{k_{x}}{\sqrt{3}}}b_{A,\mathbf{k}}^{2},e^{i\frac{k_{x}}{\sqrt{3}}}b_{B,\mathbf{k}}^{2})^{T}$.
Here $\mathbf{k}=(k_{x},k_{y})$,
\begin{equation*}
  h_{\perp,\mathbf{k}}^{0}=\left(\begin{array}{cc}
  J_{\perp} & J_{\perp} \\
  J_{\perp}e^{-i\sqrt{3}k_{x}} & J_{\perp} \end{array}
  \right)
\end{equation*}
For above bosonic Hamiltonian, it is hard to find the pseudo-time-reversal-symmetry (PTRS). But the magnon bands are still doubly degenerate when $K_{x}=0$.
This is because in this case, the bosonic Hamiltonian can be divided into two sectors and and EPA describes the coupling between the two sectors.
In two sectors, the interlayer coupling are the same with respect to a $k_{x}$-dependent phase factor. This factor coming from the interlayer glide.
For the AFM coupling, the obtained magnon bands are only dependent on the module of the coupling after Bogoliubov transformation.
So the magnon bands are doubly degenerate when $K_{x}=0$.
The EPA term couples the two sectors and will breaks the degeneracy.
In Fig.~S4 (a) and (b), we can see the EPA opens gap for the degenerate edge states and gives us second-order topological magnons, i.e., corner modes.

For the bulk with monoclinic materials, there should be 6 layers in a primitive cell to build a period along $z$-direction when considering the magnetic order.
Here to simplify our calculation, we only consider
a 2-layer period along $z$-direction. That is to say, the first and third layer are in AA-stacking. The interlayer interaction Hamiltonian
is
\begin{equation*}
  h_{\perp,\mathbf{k}}=\left(\begin{array}{cc}
  1 & 1 \\
  e^{-i\sqrt{3}k_{x}} & 1 \end{array}
  \right)\times [(J_{\perp}+\delta J_{\perp})+(J_{\perp}-\delta J_{\perp})e^{-ik_{z}}].
\end{equation*}
By considering a parallelogram pillar geometry, we can also get the hinge modes along the stacking direction, as shown in Fig.~S4 (c)
The bulk polarization along the stacking direction coming from the SSH configuration, independent on
the stacking type. The $Z_{2}\times Z$ topology is preserved and chiral hinge modes are expected to remain
in the monoclinic stacking.

\end{widetext}

\end{document}

redpajama_set_name: RedPajamaArXiv


text:
\section[1]{Introduction }
The standard cosmological model explains observations only under
the circumstance that almost all matter and energy of the Universe
are hidden in some unknown {\it dark} forms. In my model of
low-energy quantum gravity based on the idea of an existence of
the background of super-strong interacting gravitons (for more
details, see \cite{500}), a cosmological redshift is caused by
interactions of photons with gravitons. Non-forehead collisions
with gravitons lead to a very specific additional relaxation of
any photonic flux that gives a possibility of another
interpretation of supernovae 1a data - without any kinematics or
dark energy \cite{500}. I would like to summarize here the main
cosmologically essential consequences of this model. Average
energy losses of a photon with an energy $E $ on a way $dr$
through the graviton background will be equal to: $dE=-aE dr,$
where $a=H/c,$ $H$ is the Hubble constant. If we introduce a new
dimensional constant $D$, so that: $\sigma (E,\epsilon)= D \cdot E
\cdot \epsilon,$ $\sigma (E,\epsilon)$ is a cross-section of
interaction by forehead collisions of a photon with an energy $E$
and a graviton with an energy $\epsilon,$ then we can compute the
Hubble constant in this approach: $H= (1/2\pi) D \cdot \bar
\epsilon \cdot (\sigma T^{4}),$ where $\bar \epsilon$ is an
average graviton energy, and $T$ is a temperature of the
background. The constant $D$ should have the value: $D=0.795 \cdot
10^{-27}{m^{2} / eV^{2}};$ the one may be found from the Newtonian
limit of gravity. If $r$ is a geometrical distance from a source,
then we have for $r(z)$, $z$ is a redshift: $r(z)= ln (1+z)/a.$
None-forehead collisions of photons with gravitons of the
background will lead to a scatter of photons and to an additional
relaxation of a photonic flux, so that the luminosity distance
$D_{L}$ is equal in this approach to: $D_{L}=a^{-1} \ln(1+z)\cdot
(1+z)^{(1+b)/2} \equiv a^{-1}f_{1}(z),$ where $f_{1}(z)\equiv
\ln(1+z)\cdot (1+z)^{(1+b)/2}$ is the luminosity distance in units
of $c/H.$ This luminosity distance function fits supernova
observations very well for roughly $z < 0.5$. It excludes a need
of any dark energy to explain supernovae dimming.
\par
In this paper, I consider galaxy number counts/redshift and
counts/magnitude relations on a basis of this model. I assume here
that a space is flat and the Universe is not expanding.
\section[2]{The galaxy number counts-redshift relation }
Total galaxy number counts $dN(r)$ for a volume element
$dV=d\Omega r^{2}dr$ is equal to: $dN(r)=n_{g}dV=n_{g}d\Omega
r^{2}dr,$ where $n_{g}$ is a galaxy number density (it is constant
in the no-evolution scenario), $d\Omega$ is a solid angle element.
Using the function $r(z)$ of this model, we can re-write  galaxy
number counts as a function of a redshift $z$:
\begin{equation}dN(z)=n_{g}d\Omega (H/c)^{-3} {ln^{2}(1+z) \over
{1+z}}dz. \end{equation}
Let us introduce a function (see \cite{123}) $$f_{2}(z)\equiv
{(H/c)^{3}dN(z) \over n_{g}d\Omega z^{2} dz};$$ then we have for
it in this model:
\begin{equation}f_{2}(z)= {ln^{2}(1+z) \over {z^{2}(1+z)}}.\end{equation}
A graph of this function is shown in Fig. 1; the typical error bar
and data point are added here from paper \cite{72} by Loh and
Spillar. There is not a visible contradiction with observations.
{\it There is not any free  parameter in the model to fit this
curve;} it is a very rigid case.
\begin{figure}[th]
\epsfxsize=12.98cm \centerline{\epsfbox{gnc.eps}} \caption{Number
counts $f_{2}$ as a function of the redshift in this model. The
typical error bar and data point are taken from paper \cite{72} by
Loh and Spillar.}
\end{figure}
\par It is impossible to count a {\it total} galaxy number for big
redshifts so as very faint galaxies are not observable. For
objects with a fixed luminosity, it is easy to find how their
magnitude $m$ changes with a redshift. So as $dm(z)$ under a
constant luminosity is equal to: $dm(z)=5 d(lg D_{L}(z)),$ we have
for $\Delta m(z_{1},z_{2})\equiv \int_{z_{1}}^{z_{2}} dm(z):$
\begin{equation}\Delta m(z_{1},z_{2})=5 lg(f_{1}(z_{2})/f_{1}(z_{1})).
\end{equation}
This function is shown in Fig.2 for $z_{1}=0.001; 0.01; 0.1.$
\begin{figure}[th]
\epsfxsize=12.98cm \centerline{\epsfbox{magn.eps}}
\caption{Magnitude changes $\Delta m$  as a function of the
redshift difference $z_{2}-z_{1}$ in this model for $z_{1}=0.001$
(solid); 0.01 (dot); 0.1 (dash).}
\end{figure}
\par I would like to note that a very fast {\it initial} growth of
the luminosity distance with a redshift $z$ in this model might
explain the observed excess of faint blue galaxy number counts
above an expected one in the standard model (for example, see
\cite{122}). A galaxy color depends on a redshift, and a galaxy
dimming depends on the luminosity distance, because by big values
of the ratio $\Delta m(z_{1},z_{2})/(z_{2}-z_{1})$ in a region of
small redshifts and by a further much slower change of it (see
Fig.3) an observer will see many faint but blue enough galaxies in
this region (in the no-evolution scenario).
\begin{figure}[th]
\epsfxsize=12.98cm \centerline{\epsfbox{rmagn.eps}} \caption{To a
possible explanation of the excess of faint blue galaxy number
counts: $\Delta m(z_{1},z_{2})/(z_{2}-z_{1})$ vs. the redshift
difference $z_{2}-z_{1}$ in this model for $z_{1}=0.001$ (solid);
0.01 (dot); 0.1 (dash).}
\end{figure}
\section[3]{Taking into account the galaxy luminosity function}
Galaxies have different luminosities  $L,$ and we can write
$n_{g}$ as an integral: $n_{g}=\int dn_{g}(L),$ where
$dn_{g}(L)=\eta (L)dL,$ $\eta (L)$ is the galaxy luminosity
function. I shall use here the Schechter luminosity function
\cite{73}:
\begin{equation}
\eta (L)dL=\phi_{\ast}({L \over L_{\ast}})^{\alpha} exp(-{L \over
L_{\ast}}) d({L \over L_{\ast}})\end{equation} with the parameters
$\phi_{\ast},$ $L_{\ast},$ $\alpha.$\footnote{To turn aside the
problem with divergencies of this function by small $L$ for
negative values of $\alpha,$  all computations are performed here
for $z>0.001.$} So as we have by a definition of the luminosity
distance $D_{L}(z)$ that a light flux $I$ is equal to: $I={L \over
4\pi D_{L}^{2}(z)},$ and a visible magnitude $m$ of an object is
$m=-2.5 \lg I + C,$ where $C$ is a constant, then $m$ is equal to:
\begin{equation}m=-2.5 \lg I +5 \lg D_{L}(z)
+(C-4\pi).\end{equation} We can write for $L:$
\begin{equation} L=A\cdot{D_{L}^{2}(z) \over \kappa^{m}},\end{equation}
where $\kappa=10^{0.4}, \ A=const. $ For a thin layer with
$z=const$ we have:
$$ dL={\partial L \over \partial m}\cdot dm,$$ where
\begin{equation} {\partial L \over \partial m}=- m \kappa \cdot A {D_{L}^{2}(z)
\over \kappa^{m}}= -m\kappa L.\end{equation} Then
\begin{equation}
dn_{g}(m,z)=-(\phi_{\ast}\kappa)\cdot l^{\alpha}(m,z)
\exp(-l(m,z))\cdot (m\cdot l(m,z))dm,
\end{equation}
where $(-dm)$ corresponds to decreasing $m$ by growing $L$ when
$z=const,$ and $$l(m,z)\equiv {L(m,z) \over L_{\ast}}.$$
\par Let us introduce a function $f_{3}(m,z)$ with a differential
\begin{equation}
df_{3}(m,z)\equiv {dN(m,z) \over d\Omega (-dm).}
\end{equation}
We have for this differential in the model:
\begin{equation}
df_{3}(m,z)=({\phi_{\ast}\kappa \over a^3})\cdot m\cdot
l^{\alpha+1}(m,z)\cdot \exp(-l(m,z))\cdot {ln^{2}(1+z) \over
(1+z)}dz,
\end{equation}
where $a=H/c$, $H$ is the Hubble constant. An integral on $z$
gives the galaxy number counts/magnitude relation:
\begin{equation}
f_{3}(m)=({\phi_{\ast}\kappa \over a^3})\cdot m\cdot
\int_{0}^{z_{max}} l^{\alpha+1}(m,z)\cdot \exp(-l(m,z))\cdot
{ln^{2}(1+z) \over (1+z)}dz;
\end{equation}
I use here an upper limit $z_{max} =10.$ To compare this function
with observations by Yasuda et al. \cite{77}, let us choose the
normalizing factor from the condition: $f_{3}(16)=a(16), $ where
\begin{equation}
a(m)\equiv A_{\lambda}\cdot 10^{0.6(m-16)}
\end{equation}
is the function assuming "Euclidean" geometry and giving the best
fit to observations \cite{77}, $A_{\lambda}=const$ depends on the
spectral band. In this case, we have two free parameters -
$\alpha$ and $L_{\ast}$ - to fit observations, and the latter one
is connected with a constant $A_{1}\equiv {A \over a^{2}L_{\ast}}$
if $$l(m,z)=A_{1}{f_{1}^{2}(z) \over \kappa^{m}}.$$
\par If we use the magnitude scale in which $m=0$ for Vega then
$C=2.5\lg I_{Vega},$ and we get for $A_{1}$ by $H=2.14\cdot
10^{-18}\ s^{-1}$ (it is a theoretical estimate of $H$ in this
model \cite{500}):
\begin{equation}
A_{1}\simeq 5\cdot 10^{17}\cdot {L_{\odot} \over L_{\ast}},
\end{equation}
where $L_{\odot}$ is the Sun luminosity; the following values are
used: $L_{Vega}=50 L_{\odot},$ the distance to Vega $r_{Vega}=26 \
LY.$
\par Without the factor $m,$ the function $f_{3}(m)$ by $\exp(-l(m,z)\rightarrow
1$ would be close to $a(m)$ by $\alpha =-2.5.$ Matching values of
$\alpha$ shows that $f_{3}(m)$ is the closest to $a(m)$ in the
range $10<m<20$ by $\alpha =-2.43.$  The ratio ${{f_{3}(m)
-a(m)}\over a(m)}$ is shown in Fig.4 for different values of
$A_{1}$ by this value of $\alpha$. All such the curves conflow by
\begin{figure}[th]
\epsfxsize=12.98cm \centerline{\epsfbox{gnc-f3-a-2.43-A1-30.eps}}
\caption{The relative difference $(f_{3}(m)-a(m))/a(m)$ as a
function of the magnitude $m$ for $\alpha=-2.43$ by
$10^{-2}<A_{1}<10^{2}$ (solid), $A_{1}=10^{4}$ (dash),
$A_{1}=10^{5}$ (dot), $A_{1}=10^{6}$ (dadot). }
\end{figure}
$A_{1}\leq 10^2$ (or $5\cdot 10^{15}<L_{\ast}$), i.e. observations
of the galaxy number counts/magnitude relation are {\it
non-sensitive} to $A_{1}$ in this range. For fainter magnitudes
$20<m<30$, the behavior of all curves is identical: they go below
of the ratio value $1$ with the same slope. If we compare this
figure with Figs. 6,10,12 from \cite{77}, we see that the
considered model provides a no-worse fit to observations than the
function $a(m)$ if the same K-corrections are added (I think that
even a better one if one takes into account positions of
observational points in Figs. 6,10,12 from \cite{77} by $m<16$ and
$m>16$) for the range $10^2<A_{1}< 10^7$ that corresponds to
$5\cdot 10^{15}>L_{\ast}> 5\cdot 10^{10}.$
\par Observations of $N(z)$ for different magnitudes are a lot more
informative. If we define a function $f_{4}(m,z)$ as
\begin{equation}
f_{4}(m,z)\equiv ({a^3 \over \phi_{\ast}\kappa})\cdot {df_{3}(m,z)
\over dz},
\end{equation}
this function is equal in the model to:
\begin{equation}
f_{4}(m,z)= m\cdot l^{\alpha+1}(m,z)\cdot \exp(-l(m,z))\cdot
{ln^{2}(1+z) \over (1+z)}.
\end{equation}
\par Galaxy number counts in the range $m_{1}<m<m_{2}$ are
proportional to the function:
\begin{equation}
f_{5}(m_{1},m_{2})\equiv \int_{m_{1}}^{m_{2}} f_{4}(m,z)dm=
\end{equation}
$$=
 \int_{m_{1}}^{m_{2}} m\cdot l^{\alpha+1}(m,z)\cdot
\exp(-l(m,z))\cdot {ln^{2}(1+z) \over (1+z)}dm.$$
\begin{figure}[th]
\epsfxsize=10.98cm\centerline{\epsfbox{gnc-f5m,z-f4m,z10-5-10-x0.eps}}
\caption{Number counts $f_{4}(m,z)$ (dot) and $f_{5}(m_{1},m_{2})$
(solid) (logarithmic scale) as a function of the redshift by
$A_{1}=10^{5}$ for $\alpha=-2.43,\ m_{1}=10$ and  different values
of $m=m_{2}:$ 15 (blue), 20 (red), 25 (green), and 30 (mag);
$m=10$ (brown, only $f_{4}(m,z)$),.}
\end{figure}
Graphs of both $f_{4}(m,z)$ and $f_{5}(m_{1},m_{2})$ are shown in
Fig. 5 by $\alpha =-2.43,\ A_{1}=10^{5} $; they are very similar
between themselves. We see that even the observational fact that a
number of visible galaxies by $z\sim 10$ is very small allows us
to restrict a value of the parameter $A_{1}$ much stronger than
observations of $N(m)$.
\section[4]{Quasar number counts}
For quasars, we can attempt to compute the galaxy number
counts/redshift relation using  Eq. 16 with another luminosity
function $\eta{'} (l(m,z))$:
\begin{equation}
f_{5}(m_{1},m_{2})\equiv \int_{m_{1}}^{m_{2}} f_{4}{'}(m,z)dm=
 \int_{m_{1}}^{m_{2}} m\cdot l(m,z)\cdot \eta{'} (l(m,z))\cdot
{ln^{2}(1+z) \over (1+z)}dz.
\end{equation}
The following luminosity functions were probed here (see Fig. 6):
the Schechter one with $\alpha =0$, $A_{1}=10^{6.55}$ (blue); the
double power law \cite{74,76}:
\begin{equation}
\eta{'} (l(m,z)) \propto {1 \over l^{-\alpha}(m,z) +
l^{-\beta}(m,z)}
\end{equation}
with $\alpha =-3.9,\ \beta=1.6, \ A_{1}=4.5\cdot 10^{6}$ (green);
the Gaussian one:
\begin{equation}
\eta{'} (l(m,z)) \propto \exp({-(l(m,z)-1)^{2} \over 2\sigma^2})
\end{equation}
with $\sigma =0.5,\  A_{1}=4.5\cdot 10^{6}$ (brown, dot); the
combined one:
\begin{equation}
\eta{'} (l(m,z)) \propto l^{\alpha}(m,z)\cdot
\exp({-(l(m,z)-1)^{2} \over 2\sigma^2})
\end{equation}
with two sets of parameters: $\alpha =-1.45, \ \sigma =0.6, \
A_{1}=1.3\cdot 10^{6}$ (red, solid) and $\alpha =-1.4, \ \sigma
=0.7, \ A_{1}=3\cdot 10^{6}$ (red, dot). There is a couple of
\begin{figure}[th]
\epsfxsize=12.98cm
\centerline{\epsfbox{gnc-f5m,zQSOcomparison16-18.25-21.5.eps}}
\caption{QSO number counts $f_{5}(m,z)$ (arbitrary units) as a
function of the redshift for different luminosity functions:
Gaussian (brown, dot), the double power law (green), Schechter's
(blue), combined (red, solid and dot) with parameters given in the
text. The left-shifted curve of each couple corresponds to the
range $16<m<18.25,$ another one corresponds to $18.25<m<20.85.$}
\end{figure}
curves for each case: the left-shifted curve of any couple
corresponds to the range $16<m<18.25,$ another one corresponds to
$18.25<m<20.85.$ These ranges are chosen the same as in the paper
by Croom et al. \cite{74}, and you may compare this figure with
Fig. 3 in \cite{74}. We can see that the theoretical distributions
reflect only some features of the observed ones but not an entire
picture. In all these cases, a slope of an analog of
$\log(f_{3}(m))$ near $m=18$ is in the range 0.29 - 0.325, when
quasar observations give a larger slope (see Fig. 4, 21 in
\cite{74} and Fig. 13 in \cite{75}; in the latter paper, this
slope has been evaluated to be equal to about 1). We can summarize
that, as well as in the standard cosmological model, it is
impossible to fit quasar observations using some simple luminosity
function with fixed parameters.

\par In the standard model, an easy way exists to turn aside this
difficulty: one ascribes it to a quasar "evolution", then a
luminosity function (for example, the double power law
\cite{74,76}) is modified for different redshifts to take into
account this "evolution". There exist two manner to do it: one may
consider $L_{\ast}$ as a function of a redshift (pure luminosity
evolution) \cite{74} or one may assume that indices $\alpha$ and
$\beta$ of the distribution (double power law) vary with $z$
\cite{76} - in both variants, it is possible to fit observations
in some range of redshifts; of course, there are many other
descriptions of the "evolution" \cite{75}. It is strange only that
"evolutions" are not concerted: we can see exponential, quadratic
and other kinds of them - and it means that there is not any real
evolution: we deal with a pure fine art of fitting, nothing more.
In the considered model, this way is forbidden.
\par I think that it is necessary to consider some theoretical
model of a quasar activity to get a distribution of
"instantaneous" luminosities. It is known that the typical
lifetime of individual quasars is uncertain by several orders of
magnitude; a lifetime of  $4 \cdot 10^7$ years may be considered
as an average value \cite{771}. If one considers a quasar light
curve $L(t)$ (in a manner which is similar to the one by Hopkins
et al. \cite{78}) in a parametric form, it is possible to get the
luminosity function which takes into account a probability to
observe a quasar with a given luminosity. Let us consider the two
simple examples. The simplest case is a constant luminosity $L$ of
any quasar during its lifetime $\tau$. If initial moments of
quasar activity are distributed uniformly in time and may be
described by a frequency $\nu$, then a probability $P_{obs}$ to
observe a quasar will be equal to:
\begin{equation}
P_{obs}=\int_{0}^{\tau}\exp(-\nu t{'})d(\nu
t{'})=1-\exp(-\nu\tau).
\end{equation}
For $\nu\tau\ll 1$ we have $P_{obs}\simeq \nu\tau.$ If we further
assume that $\tau\propto 1/L,$ i.e. that a full emitted quasar
energy is constant, then a distribution of observable luminosities
is
\begin{equation}
\eta{'}(L)\propto \eta(L)\cdot 1/L,
\end{equation}
where $\eta(L)$ is an initial distribution of values of $L$.
\par The second example is the quasar exponential light curve:
\begin{equation}
L(t)=L_{0}\exp(-t/\tau),
\end{equation}
where $\tau$ is a lifetime, $L_{0}$ is an initial luminosity. If
$L_{0}$ has a distribution $\eta(L_{0})$, then we get:
\begin{equation}
\eta{'}(L) \propto \int
\eta(L_{0})\cdot[\exp(t_{max}/\tau(L_{0}))-1]^{-1}\cdot
(L_{0}/L)^{2-\nu\tau(L_{0})}\cdot dL_{0},
\end{equation}
where $t_{max}$ is a maximum time during which one can distinguish
a quasar from a host galaxy, and $\tau$ depends on $L_{0}$ in some
manner. We see that even in this simple toy example the dependence
on $\tau$ is not trivial.
\par In a general case, it is necessary to describe both - front
and back - slopes of a quasar light curve. Together with a total
emitted energy (or a peak luminosity), we need at least three
independent parameters; if we take into account their random
distributions, this number should be at least doubled.

\section[5]{Conclusion}
Starting from a micro level and considering interactions of
photons with single gravitons, we can find the luminosity distance
and a geometrical distance in this approach. Using only these
quantities, I compute here galaxy number counts-redshift and
galaxy number counts-magnitude relations for a case of a flat
non-expanding universe. It has been shown here that they are in a
good accordance with observations. It may be important as for
cosmology as for a theory of gravity.


text:
\section*{Introduction}
\input{intro}

\numberwithin{equation}{section}

\section[A chain of inequalities]
  {\raggedright A chain of inequalities}
\label{sect1}
\input{sect1}

\section[Splittable random processes]
  {\raggedright Splittable random processes}
\label{sect2}
\input{sect2}

\section[Logarithm of a Gaussian process]
  {\raggedright Logarithm of a Gaussian process}
\label{sect3}
\input{sect3}

\section[An example of a splittable process]
  {\raggedright An example of a splittable process}
\label{sect4}
\input{sect4}

\section[Splittable random fields]
  {\raggedright Splittable random fields}
\label{sect5}
\input{sect5}

\section[An example of a splittable field]
  {\raggedright An example of a splittable field}
\label{sect6}
\input{sect6}

\section[Random complex zeroes]
  {\raggedright Random complex zeroes}
\label{sect7}
\input{sect7}


text:
\section{Introduction}
Beamforming is a signal processing technique that utilizes the spatial information from a set of receivers to attain the signal with satisfactory quality. It has been used in different applications, wireless communications, acoustics, etc.    As a data-dependent beamformer, the classical Capon method \cite{capon1969} is well known to achieve this goal and maximizes the output signal-to-interference-plus-noise (SINR) level using the minimum variance distortionless response (MVDR) criterion.
However, in a presence of a steering vector error caused by, for example, a mismatch in a direction of arrival (DoA), the performance of the beamformer seriously degraded since the desired signal is mistreated as an interference and suffers a severe attenuation.
In order to resolve this issue, numerous approaches have been proposed for dealing with the steering vector errors.
A class of beamformers known as linear constrained minimum variance (LCMV) were proposed to mitigate this issue by imposing additional linear constraints that increase the beamwidth of the mainlobe \cite{Vural1977,Takao1976}.
Another popular class of beamformers, called diagonal loading (DL), addressed the issue by adding a scalar matrix to the covariance matrix \cite{bd1998,li2003}.
DL methods are known to provide robustness not only to the DoA mismatch issue but also to general types of steering vector errors.
One well-known drawback of the DL methods is that a universally accepted way of choosing the diagonal loading factor is still lacking.
Nevertheless, in \cite{vincent2004}, an optimal DL factor that maximizes SINR was derived with some approximations.
Further, in \cite{ma2003,du2010}, automatic determination of DL factor values was presented.
In order to enhance the performance, many extensions to DL \cite{sq1999,vorobyov2003, Lorenz2005,ShahbazPanahi2003,chen2007} were proposed with an idea to force the magnitude response of all steering vectors in a sphere set centered at the desired steering vector to exceed or equal to unity while minimizing the output variance.
In \cite{Feldman1994}, a method is proposed by projecting the steering vector onto the signal-plus-interference subspace. More recently, methods that based on complex quaternion processes \cite{tao2014}, based on the artificial immune system \cite{kiong2014}, based on interference-plus-noise covariance matrix reconstruction \cite{huang2015}, based on an interference subspace \cite{mao2015} were also proposed.

In this study, we consider only the steering vector error caused by DoA mismatches and propose a new extension to DL method called subspace-constrained diagonal loading (SSC-DL).
The difference between SSC-DL and conventional DL is that the search space of weight vectors is constrained in a subspace determined by the signal covariance matrix and the possible range of the desired signal's DoA.
An important feature of the proposed method over other DL methods is that the additional constraints are added to ensure an interference suppression capability of the beamformer instead of working only on broadening the beamwidth.
Such constraints turn out to be helpful in improving SINR performance.
The simulation results show the effectiveness of the proposed SSC-DL approach in suppressing interference signals.

The rest of this paper is organized as follows.
Section II introduces the MVDR beamforming system and the DoA mismatch issue.
Section III presents the proposed SSC-DL algorithm.
Section VI presents the derivation of the diagonal loading factor used in the proposed SSC-DL algorithm.
Section \ref{sec:sim} presents the simulation setup and results.
Finally, conclusions are provided in Section \ref{sec:conclusion}.
Boldface lowercase letters represent column vectors; boldface uppercase letters represent matrices. {Superscripts $^T$, $^H$, $^\perp$, and $E[\cdot]$ represent matrix transpose, Hermitian-transpose operations, and subspace orthogonal complement, and expectation operation, respectively.}

\section{MVDR Beamformers}
In this section, we first introduce the beamforming system, MVDR beamformer, and the DL method. Next, we describe the DoA mismatch issue.

\subsection{The Beamforming System}
\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth,clip]{system_model.jpg}}
\caption{system model for robust beamforming}\label{fig:system}
\end{figure}

With figure \ref{fig:system}, consider a beamforming system comprising $N$ sensors in a uniform linear array (ULA) with an inter-element spacing $d$.
The estimated DoA of the desired signal is $\theta_0$, and the actual DoA of the desired signal is $\theta_d$.
When $\theta_0$ is not equal to $\theta_d$, the difference $\theta_0-\theta_d$ is called the DoA mismatch.
We denote $\mathbf{a}(\theta)\in\mathbb{C}^N$ as the steering vector toward direction $\theta$ \cite{bd1998}: \begin{align}\label{eq2}
\mathbf{a}(\theta)=[1,e^{j(2\pi/\lambda)d\sin\theta},\cdots,e^{j(N-1)(2\pi/\lambda)d\sin\theta}]^{T}
\end{align}
where $\lambda$ is the wavelength.
Assume there are $J$ interference signals $x_j(t), j=1,..., J$, and the $j$th interference signal is arriving from angle $\theta_{Ij}$, {where $t\in{\cal R}$ denotes time.}
Then, the received vector signal  {at time $t$}, $\mathbf{y}(t)=[y_1(t), y_2(t),\cdots,y_N(t)]^T$ can be expressed as
\begin{equation}
\begin{aligned}\label{eq:yt}
\mathbf{y}(t)&=x_0(t)\mathbf{a}(\theta_d)+\sum_{j=1}^J x_{j}(t)\mathbf{a}(\theta_{Ij})+\mathbf{n}(t)\\
&=x_0(t)\mathbf{a}(\theta_d)+{\mathbf i}(t)+{\mathbf n}(t)
\end{aligned}
\end{equation}
where $x_0(t)$ is the desired signal and {${\bf n}(t)$ denotes additive white Gaussian noise whose entries are i.i.d. complex Gaussian random variables $\sim {\cal CN}(0, \sigma_n^2)$. The signals $x_0(t)$ and $x_{j}(t), j=1,..., J$ are assumed to be statistically independent zero-mean wide-sense stationary processes, with $E[|x_0(t)|^2] = \sigma^2_{x}$ and $E[|x_{j}(t)|^2] = \sigma^2_{I_{j}}$.}
The desired signal $x_0(t)$ and the interference signals $x_j(t)$ are assumed to be uncorrelated.

The beamforming technique uses a weight vector $\mathbf{w}=[w_0, w_1,\cdots,w_{N-1}]^T$ to process the received signal ${\bf y}(t)$ in an attempt to suppress the interference component present in ${\bf y}(t)$ while preserving the desired signal $x_0(t)$.
The performance of the processed signal $  s(t) = {\bf w}^H {\bf y}(t)$ can be evaluated by the signal-to-interference-plus-noise ratio (SINR), defined as
\begin{align}\label{eq:SINR}
SINR=\frac{P\left | \mathbf{w}^H {\bf a}(\theta_d) \right |^2}{\mathbf{w}^H\mathbf{R}_{I+n}\mathbf{w}}
\end{align}
where $P$ = $E[|x_0(t)|^2]$ is the desired signal power and $\mathbf{R}_{I+n}$ = $\mathbf{R}_{I}$ + $\mathbf{R}_{n}$ is the sum of the interference-covariance matrix {$\mathbf{R}_{I}$ = $E[\mathbf{i}(t)\mathbf{i}(t)^H]$, with ${\bf i}(t)$ defined as ${\bf i}(t)=\sum_{j=1}^J x_{j}(t)\mathbf{a}(\theta_{Ij})$}, and noise-covariance matrix $\mathbf{R}_{n}$ = $E[\mathbf{n}(t)\mathbf{n}(t)^H] = \sigma_n^2{\bf I}_n$.

\subsection{MVDR Beamformer}
The MVDR beamformer estimates the weight vector $\mathbf{w}$ by solving the following optimization problem:
\begin{align}\label{eq:MVDR-RI+n}
&\min_{\mathbf{w}}\mathbf{w}^{H}\mathbf{R}_{I+n}\mathbf{w}  ~~~ s.t. ~ |\mathbf{w}^{H}\mathbf{a}(\theta_d)|=1 .
\end{align}
Generally, $\mathbf{R}_{I+n}$ is difficult to be precisely estimated, and thus the classical MVDR beamformer uses signal covariance matrix $\mathbf{R}_{y}=E[\mathbf{y}(t)\mathbf{y}(t)^H]$ to replace $\mathbf{R}_{I+n}$. {When $\theta_d$ is perfectly known and $|\mathbf{w}^{H}\mathbf{a}(\theta_d)|=1$ is satisfied}, it is shown \cite{capon1969} that (\ref{eq:MVDR-RI+n}) is equivalent to:
\begin{align}\label{eq:MVDR}
	&\min_{\mathbf{w}}\mathbf{w}^{H}{\mathbf{R}}_{y}\mathbf{w}  ~~~ s.t. ~ |\mathbf{w}^{H}\mathbf{a}(\theta_d)|=1.
\end{align}
{The equivalence of (\ref{eq:MVDR-RI+n}) and (\ref{eq:MVDR}) becomes clear by observing that
\begin{align}
	\mathbf{w}^{H}{\mathbf{R}}_{y}\mathbf{w} &= \mathbf{w}^{H}(P{\bf a}(\theta_d){\bf a}^H(\theta_d)+{\mathbf{R}}_{I+n})\mathbf{w} \nonumber\\
	&=P + \mathbf{w}^{H}{\mathbf{R}}_{I+n}\mathbf{w} .
\end{align}}
The optimal solution for (\ref{eq:MVDR}) can be expressed in a closed form:
\begin{align}\label{eq:MVDR-closedForm}
\mathbf{w}=\frac{{\mathbf{R}}_{y}^{-1}\mathbf{a}(\theta_d)}{\mathbf{a}^H(\theta_d){\mathbf{R}}_{y}^{-1}\mathbf{a}(\theta_d)}.
\end{align}
In practice, signal covariance matrix ${\bf R}_{ y}$ is obtained by averaging a number of snapshots $K$:
\begin{align}\label{eq:hatRy}
{\hat{\mathbf{R}}}_{y} = \frac{1}{K}\sum_{k=1}^{K}\mathbf{y}(k)\mathbf{y}^{H}(k).
\end{align}

\subsection{The Diagonal Loading Method}
The MVDR beamformer is known to have poor performance if $\theta_d$ is not known accurately.
When a DoA mismatch occurs, the desired signal may be mistreated as an interference and will be suppressed, thus decreasing the output SINR.
The DL method \cite{bd1998, li2003} is known to be capable of mitigating the performance degradation caused by the DoA mismatch issue by including a scalar matrix in Eq. ($\ref{eq:MVDR}$):
\begin{align}\label{eq:DL}
{\bf w}_{DL}=\arg\min_{{\bf w}\in\mathbb{C}^N} {\bf w}^H ({\bf R}_y + \gamma{\bf I}){\bf w} ~~s.t.~ |{\bf w}^H {\bf a}(\theta_0)|=1
\end{align}
where $\gamma$ is a real number.

The DL method is successful in avoiding the beamformer to misrecognize the desired signal as an interference;
it is, however, limited in the capability of suppressing interference signals, leaving room for improvements in SINR performance.

\subsection{Problem Statement}
We consider a beamforming system with a DoA mismatch: the estimated DoA of the desired signal is $\theta_0$, which is unequal to the true DoA, $\theta_d$. Moreover, we assume the receiver knows a lower bound $\theta_1$ and an upper bound $\theta_2$ of the desired DoA, i.e., $\theta_1 < \theta_d < \theta_2$. {The DoAs of interference signals, $\theta_{Ij} $, are assumed not to be within the range $[\theta_1, \theta_2]$}. {We assume that the angles $\theta_d, \theta_{Ij},j=1,...,J$ are deterministic yet in general unknown to the receiver.} The goal is to use $\hat{\bf R}_y$ defined in (\ref{eq:hatRy}) and find a weight vector ${\bf w}$ that achieves a good SINR performance defined in (\ref{eq:SINR}). Without loss of generality, we assume all angles $\theta_d$, $\theta_0$, $\theta_1, \theta_2$, and $\theta_{Ij}$ are within the range $(-\pi/2,\pi/2)$.

\section{The Proposed SSC-DL Beamformer}
Although the MVDR with DL is capable of mitigating the unwanted suppression of the desired signal, there is still room to further improve the SINR performance in the presence of the DoA mismatch.
In this section, we first describe the implementation steps of the proposed SSC-DL beamformer followed by an explanation of its advantages on interference suppression.

\subsection{Algorithm Description}
\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth,height=5cm]{SSCDL_flow_chart.jpg}}
\caption{Flow Chart of SSC-DL}
\label{flsscdl}
\end{figure}
{Figure \ref{flsscdl} shows the flow chart of the proposed SSC-DL algorithm, which can be divided into three steps.} In the first step, SSC-DL chooses an integer $M<N$ and defines a set of $M$ angles neighboring to the estimated DoA, $\theta_0$, within the range $[\theta_1, \theta_2]$:
\begin{align}\label{eq:phi_m}
\phi_m = \sin^{-1}\left(\sin{\theta_1} + m \cdot \frac{\sin{\theta_2} - \sin{\theta_1}}{M-1}\right),
\end{align}
$m=0,1,...M-1$. Note that $\theta_1=\phi_0<\phi_1<\cdots<\phi_{M-1}=\theta_2$. Then, for the $m$th angle $\phi_m$, we compute its corresponding weight vector according to (\ref{eq:MVDR-closedForm}):
\begin{align}\label{eq:w_m}
{\bf w}_m = {  \frac{{\bf R}_y^{-1}{\bf a}(\phi_m)}{{\bf a}^H(\phi_m){\bf R}_y^{-1}{\bf a}(\phi_m)}},
\end{align}
which is the optimal weight vector in terms of SINR when the true DoA happens to be $\theta_d=\phi_m$.
In the second step, SSC-DL collects these $M$ vectors
and constructs an $M$-dimensional subspace of $\mathbb{C}^N$
\begin{equation}
{\cal W} = \left\{ {\bf W}_{M}{\bf z} | {\bf z}\in\mathbb{C}^M \right\}
\label{eq:subspaceW}
\end{equation}
where ${\bf W}_{M} = \left[\begin{array}{cccc}{\bf w}_1 & {\bf w}_2 & \cdots & {\bf w}_M \end{array}\right]$ is of size $N\times M.$
The third step is to compute the weight vector as:
\begin{align}\label{eq:VS}
{\bf \hat w}_{M} = \arg\min_{{\bf w}\in {\cal W}} {\bf w}^H ({\bf R}_y + \gamma{\bf I}){\bf w} ~~~s.t. ~|{\bf w}^H {\bf a}(\theta_0)|=1.
\end{align}
Here, the value $\gamma$ is chosen as
\begin{equation} \gamma = -(\sigma_{n}^2 + P\cdot N),\label{eq:gamma}
\end{equation}
based on a derivation that will be later elaborated in Section \ref{sec:derivation}.
Using a technique of Lagrange multipliers, the optimization problem defined in (\ref{eq:VS}) can be shown to have a closed-form solution
\begin{align}\label{ssc-dl}
{\mathbf {\hat w}}_{M}=\frac{{\bf W}_{M}({\bf W}_{M}^{H}(\mathbf{R}_{y}+\gamma\mathbf{I}){\bf W}_{M})^{-1}{\bf W}_{M}^{H}\mathbf{a}(\theta_0)}{\mathbf{a}^{H}(\theta_0){\bf W}_{M}({\bf W}_{M}^{H}(\mathbf{R}_{y}+\gamma\mathbf{I}){\bf W}_{M})^{-1}{\bf W}_{M}^{H}\mathbf{a}(\theta_0)}.
\end{align}

\subsection{Choices of the subspace dimension $M$}\label{subsec:M}
The choice of the dimension of the subspace ${\cal W}$, $M$, could be any integer between $1$ and $N$.
When $M=1$, the SSC-DL reduces to the  MVDR beamformer.
On the other hand, when $M$ is chosen as $N$, the SSC-DL degenerates to the case of conventional DL method since the subspace ${\cal W}$ equals to the whole space $\mathbb{C}^M$ and there are effectively no constraints to the beamformer coefficients ${\bf w}$.

While we are not yet able to give a comprehensive criterion to an optimal choice of $M$ in this work, we provide here some rules of thumbs of the selection of the parameter.
First of all, $M$ should be chosen to be smaller than or equal to $N-J$ since the subspace ${\cal W}$ should be designed to be ``almost orthogonal'' to the interference subspace which has a dimension of $J$.
On the other hand, $M$ should be chosen sufficiently large to cover all steering vectors covering all possible DoAs of the desired signal defined as the range $[\theta_1, \theta_2]$.
In Section \ref{subsec:SINRvsM}, the effect on the system performance of the parameter $M$ will be studied through numerical results.

\subsection{Interference Suppression}
\label{subsec:interf_suppr}
We give a qualitative explanation to why the proposed SSC-DL method possesses a capability of suppressing interference signals superior to the conventional DL method.
The subspace ${\cal W}$ as defined in (\ref{eq:subspaceW}) can be interpreted as space that is ``almost orthogonal'' to all interfering steering vectors ${\bf a}(\theta_{Ij}), j=1,..., J$. It is also a space that the steering vector of the desired signal, ${\bf a}(\theta_d)$, with any possible angle $ \theta_d\in [\theta_1,\theta_2]$, would have most of its energy within.
Let ${\bf Q}_s$ and ${\bf Q}_n$ be matrices whose columns, respectively, form an orthonormal basis for the subspaces ${\cal W}$ and ${\cal W}^\perp$.
Then, the steering vector of an interfering signal will have most of its energy in ${\cal W}^\perp$:
\begin{equation} || {\bf Q}_s^H {\bf a}(\theta_{Ij})|| \ll ||{\bf Q}_n^H  {\bf a}(\theta_{Ij})|| \approx ||{\bf a}(\theta_{Ij})||,~ j=1, 2, ..., J.
\nonumber \label{eq:QsllQn}
\end{equation}
Denote $\lambda_j$ as the ratio $\lambda_j = || {\bf Q}_s^H {\bf a}(\theta_{Ij})|| / ||{\bf Q}_n^H  {\bf a}(\theta_{Ij})||$.
Then we have $\lambda_j \ll 1$.
Further, noting that ${\bf \hat w} = {\bf Q}_s{\bf z}'$ for some ${\bf z}'$ satisfying $||{\bf \hat w}|| = ||{\bf z}'||$, we have
\begin{eqnarray}
 |{\bf a}^H(\theta_{Ij}){\bf w}| &=& |{\bf a}^H(\theta_{Ij}) {\bf Q}_s {\bf z'}|  \leq  || {\bf Q}_s^H {\bf a}(\theta_{Ij})|| \cdot ||{\bf z'}|| \nonumber\\
 &=&  \lambda_j|| {\bf Q}_n^H {\bf a}(\theta_{Ij})||\cdot  || {\bf w}||\nonumber\\
&\ll&   ||{\bf a}(\theta_{Ij})|| \cdot || {\bf w}||.\nonumber
\end{eqnarray}
This means, the beamformer ${\bf \hat w}$ produced by SSC-DL guarantees a suppression level on the signal component ${\bf a}(\theta_{Ij})$ with a response below $\lambda_j$,
whose value is roughly inversely proportional to the root-mean strength of the interference ${\bf a}(\theta_{Ij})$ present in the autocorrelation matrix ${\bf R}_y$.

On the other hand, a beamformer obtained by the conventional DL method, ${\bf w}_{DL}$, does not have such a guarantee.
Adjusting the diagonal loading factor $\gamma$ in the DL method has the effect of reducing attenuation of any signal components from directions other than $\theta_0$, and is, therefore, helpful in relieving undesired suppression of the target signal from angle $\theta_d\neq\theta_0$.
However, it also reduces suppression on any interference signal from angle $\theta_{Ij}$.
By merely adjusting $\gamma$, the DL method is unable to distinguish interference signals from the desired signal and apply different suppression levels on them.
As the DL objective function (\ref{eq:DL}) does not explicitly prohibit the beamformer to contain signal components in ${\cal W}^\perp$, ${\bf w}_{DL}$ may comprise components from both ${\cal W}$ and ${\cal W}^\perp$.
Assume ${\bf w}_{DL}$ is decomposed as
$ {\bf w}_{DL} = {\bf Q}_s {\bf z}_s + {\bf Q}_n {\bf z}_n$.
Then, almost certainly ${\bf w}_{DL}$ has a nonzero (though small) component in ${\cal W}^\perp$ (i.e., $||{\bf z}_n||>0$).
Then $ | {\bf a}^H(\theta_{Ij}) {\bf \hat w}| < |{\bf a}^H(\theta_{Ij}){\bf w}_{DL}|$ is very likely to happen, as long as $||{\bf z}_n||/ ||{\bf z}_s|| \gg \lambda_j$, which is usually true, especially when ${\bf a}(\theta_{Ij})$ represents a strong interference source. { In Section \ref{secsec:ProjectionRatio}, we will provide some numerical results to demonstrate that beamformers obtained by conventional DL methods have a larger proportion of energy in the interference domain than those obtained by SSC-DL}.
The observation prevents DL methods from having a good overall SINR performance.

\section{Derivation of DL factor for SSC-DL}
\label{sec:derivation}
In this section, we present the formulation to determine the DL factor $\gamma$ for the SSC-DL method. The derivation basically follows the flow used in \cite{vincent2004}, with a major difference that the matrix ${\bf W}_M$ is involved here due to the subspace constraint. For notational simplicity, we use ${\bf a} = {\bf a}(\theta_{d})$, $\bar{\bf a} = {\bf a }(\theta_{0})$, and $\sigma_{n}^2$ is noise variance. Our goal is to determine $\gamma$ that maximizes the output SINR, i.e.,
\begin{align}\label{maxsinr}
    \gamma_{opt} = \arg\max_{\gamma}\frac{P\left | \hat{{\mathbf w}}_{M}^{H} {\bf a}(\theta_{d}) \right |^2}{\hat{{\mathbf w}}_{M}^{H}{\bf R}_{I+n}\hat{{\mathbf w}}_{M}}.
\end{align}
Substituting (\ref{ssc-dl}) into the objective function in (\ref{maxsinr}), and with some efforts similar to that in \cite{vincent2004}, it can be shown that:
\begin{align}\label{SINR_newform}
&\frac{P\left | \hat{{\mathbf w}}_{M}^{H} {\bf a}(\theta_{d}) \right |^2}{\hat{{\mathbf w}}_{M}^{H}{\bf R}_{I+n}\hat{{\mathbf w}}_{M}} = \frac{1}{P({\bf a}-\alpha\bar{{\bf a}})^H{\bf W}_{M}{\bf Z}{\bf W}_{M}^H({\bf a}-\alpha\bar{{\bf a}})}
\end{align}
where the value $\alpha$ and the matrix ${\bf Z}$ are defined as
\begin{equation}\label{alpha}
 \alpha = \frac{1+P{\bf a}^H{\bf W}_{M}({\bf W}_{M}^H({\bf R}_{I+n}+\gamma{\bf I}){\bf W}_{M})^{-1}{\bf W}_{M}^H{\bf a}}{P{\bf a}^H{\bf W}_{M}({\bf W}_{M}^H({\bf R}_{I+n}+\gamma{\bf I}){\bf W}_{M})^{-1}{\bf W}_{M}^H\bar{{\bf a}}}
\end{equation}

\begin{align}\label{Z}
\mbox{and }~~~{\bf Z} &=& ({\bf W}_{M}^H({\bf R}_{I+n}+\gamma{\bf I}){\bf W}_{M})^{-1}{\bf W}_{M}^H{\bf R}_{I+n}{\bf W}_{M}\nonumber\\
&&\cdot({\bf W}_{M}^H({\bf R}_{I+n}+\gamma{\bf I}){\bf W}_{M})^{-1}.
\end{align}
Define the function $f(\gamma)$ as the denominator of the right-hand side of Eq. (\ref{SINR_newform}):
\begin{equation}\label{f}
f(\gamma) = ({\bf a}-\alpha\bar{{\bf a}})^H{\bf W}_{M}{\bf Z}{\bf W}_{M}^H({\bf a}-\alpha\bar{{\bf a}}).
\end{equation}
Then minimizing $f(\gamma)$ is equivalent to the problem in (\ref{maxsinr}).
To simplify the expressions of $\alpha$ and ${\bf Z}$, we assume that ${\bf W}_M$ is ``almost orthogonal'' to any interfering signal ${\bf a}(\theta_j)$ (see Section \ref{subsec:interf_suppr}) and recognize that a common term in $\alpha$ and ${\bf Z}$ can be approximated as
\begin{eqnarray}
&& {\bf W}_M({\bf R}_{I+n}+\gamma{\bf I}){\bf W}_M \nonumber\\
&=& {\bf W}_M\left[\sum_{j=1}^J x_j(t){\bf a}(\theta_j){\bf a}^H(\theta_j)+(\sigma_n^2+\gamma){\bf I}\right]{\bf W}_M \nonumber\\
&\approx& (\sigma_n^2 + \gamma){\bf W}_M^H{\bf W}_M = (\sigma_n^2 + \gamma){\bf V}{\boldsymbol\Sigma}^2_M{\bf V}^H  \label{near_orthogonal}
\end{eqnarray}
where the matrix ${\bf W}_M$ is assumed to have the form of singular value decomposition
$ {\bf W}_M = {\bf Q}_s{\boldsymbol\Sigma}_M{\bf V}^H$.
Note that ${\bf Q}_s$ forms an orthonormal basis of the constraint subspace in SSC-DL.
Then, using (\ref{Z})(\ref{near_orthogonal}), we have
\begin{align}\label{WZW}
{\bf W}_M{\bf Z}{\bf W}_M^H \approx \frac{\sigma_n^2}{(\sigma_n^2 + \gamma)^2}{\bf Q}_s{\bf Q}_s^H.
\end{align}
Similarly, Eq. (\ref{alpha}) can be simplified to
\begin{align}\label{n_alpha}
&\alpha \approx  \underbrace{\frac{\sigma_n^2 + P {\bf a}_n^H {\bf a}_n}{P {\bf a}_n^H \bar{\bf a}_n}}_{\triangleq\alpha_0} +
{\gamma}\cdot\underbrace{\frac{1}{P {\bf a}_n^H \bar{{\bf a}}_n}}_{\triangleq\alpha'_0}
\end{align}
where ${\bf a}_{n} = {\bf Q}_s^H{\bf a}$ and $ \bar{{\bf a}}_{n} = {\bf Q}_s^H\bar{{\bf a}} $ represent projections of desired and estimated steering vectors ${\bf a}$ and $\bar{\bf a}$ onto the constraining subspace. Then, using (\ref{WZW}), (\ref{n_alpha}), the function (\ref{f}) can be rewritten as
\begin{align}\label{fgamma}
f(\gamma) = \frac{\sigma_n^2}{(\sigma_n^2 + \gamma)^2} \|{\bf b}+\gamma{{\bf b}}'\|^2,
\end{align}
where
\begin{align}\label{algebra}
&{\bf b} = {\bf Q}_s^H ({\bf a}-\alpha_0\bar{{\bf a}}) = {\bf a}_{n} - \alpha_{0}\bar{{\bf a}}_{n} \\
&{\bf{b}}' = -{\bf Q}_s^H{\alpha_0}'\bar{{\bf a}} = -{\alpha_0}'\bar{{\bf a}}_{n} \nonumber
\end{align}
Then, by taking derivative on $f$ in (\ref{fgamma}) with respective to $\gamma$, the optimal solution of $\gamma$ can be found to be
\begin{equation}\label{opt_gamma}
\gamma_{opt} = \frac{\|{\bf b}\|^2 - \sigma_n^2\mathbb{R}\left [ {\bf b}^H{{\bf b}}'\right ]}{\sigma_n^2\|{\bf{b}}'\|^2-\mathbb{R}\left [ {\bf b}^H{\bf{b}}'\right ]}
= - (P\|{\bf a}_n\|^2 + \sigma_n^2).
\end{equation}
Note that the norm of ${\bf a}_n={\bf Q}_s^H{\bf a}$ would be very close to that of ${\bf a}$ as long as the  steering vector of the desired signal has most of its energy in the constraining subspace, leading to
\begin{equation}
\|{\bf a}_{n}\|^2 \cong \|{\bf a}\|^2 = N.\label{ana}
\end{equation}
So (\ref{opt_gamma}) becomes
\begin{equation}\label{hatgamma}
\hat{\gamma} = -(\sigma_n^2 + P\|{\bf a}\|^2) =  -(\sigma_n^2 + P\cdot N).
\end{equation}
This result happens to be identical to that given in \cite{vincent2004}. However, the validity of this optimal solution is subject to two assumptions given in (\ref{near_orthogonal}) and (\ref{ana}). As in SSC-DL, these two assumptions are satisfied, the choice of $\gamma$ presented in (\ref{eq:gamma}) is justified.

\section{Simulation}\label{sec:sim}
In this section, we report the simulation results of the proposed SSC-DL algorithm.
For comparison purpose, we conducted experiments on several related beamforming algorithms, including (a) MVDR\cite{capon1969}, (b) MVDR with DL \cite{bd1998, li2003, du2010}, (c) a quadratically constrained method with DoA beamwidth widening in \cite{chen2007}, and (d) an autocorrelation matrix reconstruction method in \cite{huang2015}.
These four methods are labeled as ``MVDR'', ``MVDR-DL'', ``[Chen2007],'' and ``[Huang2015],'' respectively, in the following discussions.
The optimal result of MVDR method with perfect $\mathbf{R}_{y}$ and without any DoA error is also reported and denoted as ``Optimal,'' which serves as the theoretical limit for all methods.

\subsection{Simulation Setup}
\label{subsec:simsetup}
{In our simulation}, the uniform linear array is employed with $N=10$ omni-directional sensors whose inter-element distance is half the wavelength: $\lambda/2$.
Except for ``Optimal'', all methods use $K = 100$ snapshots to construct the $\hat{\bf{R}}_y$ matrix defined in (\ref{eq:hatRy}).
The actual DoA of the desired signal is set to $\theta_d=0^{\circ}$ and the estimated DoA is $\theta_0=2.5^{\circ}$.
Two interference sources are radiating through the angles $\theta_{I1} = -20^\circ$ and $\theta_{I2} = 30^\circ$ (i.e., $J=2$). Then, ($\ref{eq:yt}$) is rewritten as
\begin{equation*}\label{eq20}
\mathbf{y}(t)=x_{0}(t)\mathbf{a}(\theta_{d})+x_{1}(t)\mathbf{a}(\theta_{I1})+x_{2}(t)\mathbf{a}(\theta_{I2})+\mathbf{n}(t).
\end{equation*}
{The $x_{0}(t)$, $x_{1}(t)$, $x_{2}(t)$, and $\mathbf{n}(t)$ are zero mean wide-sense independent stationary random processes and satisfy}:
\begin{eqnarray*}
E[|x_{0}(t)|^{2}]=\sigma_{x}^{2}= P, &&
E[|x_{1}(t)|^{2}]=\sigma_{I_{1}}^{2}=10^4(40\mbox{dB}) \\
E[\mathbf{n}(t)\mathbf{n}^{H}(t)]=\sigma_{n}^{2}\mathbf{I}_{N}, &&
E[|x_{2}(t)|^{2}]=\sigma_{I_{2}}^{2}=10^2(20\mbox{dB})
\end{eqnarray*}
where $\sigma_{n}^{2}$ = 1. The perfect covariance matrix is expressed as
\begin{equation}
\begin{aligned}\label{eq21}
\mathbf{R}_{y}&=P\mathbf{a}(\theta_{d})\mathbf{a}^{H}(\theta_{d})+\sigma_{I_{1}}^{2}\mathbf{a}(\theta_{I1})\mathbf{a}^{H}(\theta_{I1})\\
&+\sigma_{I_{2}}^{2}\mathbf{a}(\theta_{I2})\mathbf{a}^{H}(\theta_{I2})+\mathbf{I}_{N}.
\end{aligned}
\end{equation}
All beamformers know that the estimated DoA $\theta_0$ has an error up to $\pm 4^\circ$ (i.e., $\theta_1 = -1.5^\circ$ and $\theta_2 = 6.5^\circ$) \cite{chen2007}.
We choose $M=5$ for the proposed SSC-DL method.
For both the MVDR-DL and the proposed SSC-DL methods, the DL factor is chosen as $\gamma = -(\sigma_n^2 + P \cdot N)$ \cite{vincent2004}.
For [Chen2007], we follow their Algorithm 2 with the initial DL factor set as $\gamma=1$, and with step angles $\zeta_1=1^\circ, \zeta_2 = 2.5^\circ, $ and $ \zeta_3 = 4.5^\circ$ \cite{chen2007}.
For [Huang2015], we follow the implementation steps reported in \cite{huang2015} by the estimated DoA of interference signals are $\bar{\theta}_{I1} = -17.5^{\circ}$ and $\bar{\theta}_{I2} = 32.5^{\circ}$, and thus $\Theta_{int}$ = $[\bar{\theta}_{I1}-4,\bar{\theta}_{I1}+4]$ $\cup$ $[\bar{\theta}_{I2}-4,\bar{\theta}_{I2}+4]$; the spherical uncertainty constant $\epsilon = 0$, the number of sampling points in $\Theta_{int}$ is $Q$ = 40, and the number of sampling points on the surface of the spherical uncertainty set is $L$=0 \cite{huang2015}.
The results of above beamforming methods, to be reported below, are the average of 200 Monte Carlo simulations with $x_{0}(t)$, $x_{1}(t)$ and $x_{2}(t)$ having the standard normal distribution.

\subsection{Magnitude Response Versus Incoming Angle}
\label{subsec:caponspectrum}
\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth,height=7cm]{magnitude_response_vs_angle.eps}}
\caption{Magnitude response $|\mathbf{w}^H\mathbf{a}(\theta)| $ versus angle $\theta$}
\label{fig4}
\end{figure}

Figure \ref{fig4}(A) shows the beamforming patterns of all methods by plotting the magnitude response versus azimuth angles in degrees.
The SNR $P/\sigma_n^2$ is set to 10 dB in this plot.
Since the desired signal DoA is at $\theta_d=0^{\circ}$ (red vertical line), and the two interference sources are at $\theta_{I1} = -20^\circ$ and $\theta_{I2} = 30^\circ$(black vertical lines), a successful beamforming technique should give high magnitude response around $0^{\circ}$ and low magnitude response around $-20^\circ$ and $30^\circ$.
Figures \ref{fig4}(B)(C)(D) display zoom-in looks of Figure \ref{fig4}(A) at angles around $\theta_d$, $\theta_{I1}$, and $\theta_{I2}$, respectively.
From Figure \ref{fig4}(B), we note that due to the DoA mismatch, MVDR gives low magnitude response at $\theta_d=0^{\circ}$ while all  other methods are robust to such a mismatch and generate close-to-unity magnitude responses at the desired DoA.
From Figures \ref{fig4}(C)(D), SSC-DL shows effective interference signal suppression performance, especially at the angle of $-20^\circ$ where a strong interference resides.
At the angle of $30^\circ$, SSC-DL does not have a lower response than [Chen2007] and [Huang2015].
But the interferer at $30^\circ$ is 20dB weaker than that at $-20^\circ$, so SSC-DL turns out to still have the largest SINR among all methods, as will be validated later in Section \ref{sub:SINRvsSNR}.
Another interesting observation is the interference signal magnitude responses  of SSC-DL are very close to that of Optimal.

Further, we compare the main lobes of beamformers produced by SSC-DL and [Chen2007].
Both methods are given the information that the desired angle $\theta_d$ is within the range $[\theta_1, \theta_2] = [-1.5^\circ, 6.5^\circ]$.
The method [Chen2007] uses a strategy to force magnitude responses for all angles in $[-1.5^\circ, 6.5^\circ]$ to be greater than or equal to unity, as can be seen from Figure \ref{fig4}(B).
The SSC-DL, on the other hand, has a narrower beamwidth, with large magnitude responses only at around the desired signal's true DoA.
This suggests setting explicit constraints on maintaining magnitude responses at the range of all possible angles is not necessary.

\subsection{Output SINR versus SNR}

\label{sub:SINRvsSNR}
\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{7}}
\caption{Output SINR versus SNR  }
\label{fig10}
\end{figure}

Figure \ref{fig10} shows the correlation of input SNR ($P/\sigma_n^2$) to the output SINR.
The plot is generated by fixing $\sigma_n^2, \sigma_{I1}^2$, and $\sigma_{I2}^2$ and by adjusting $P$ to reflect different input SNR levels.
For this set of simulation, given an input SNR, a higher output SINR indicates a better performance.
From Figure \ref{fig10}, we first observe that MVDR-DL performs poorly when the power of the desired signal is close to the power of some interference (e.g., SNR$\approx 20$ dB); this is caused by the fact that the interference signal is not effectively suppressed, and thus the output SINR is degraded.
It is clear that SSC-DL successfully overcomes such a limitation by using the subspace technique.
Moreover, we note that SSC-DL not only outperforms all other algorithms across different input SNR levels between 0 to 30 dB but also achieves a performance very close to the Optimal benchmark.

\subsection{Output SINR versus snapshot number}
\label{subsec:SINRvsK}
\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{10}}
\caption{Output SINR versus snapshot number ($K$)}
\label{fig:SINRvsK}
\end{figure}
In Figure \ref{fig:SINRvsK}, we investigate the correlation of the number of snapshots $K$ used in producing $\hat{\bf R}_y$ defined in (\ref{eq:hatRy}) and the output SINR performance of all beamforming methods.
In addition, a new curve labeled ``MVDR, no mismatch'' is added to show the performance of method \cite{capon1969} without DoA mismatch (i.e., resetting $\theta_0=0^\circ$).
The SNR is set to 10dB in this plot.
Inheriting the good property of MVDR-DL against the finite sample effect, SSC-DL also has a good performance with a limited number of snapshots.
From Figure \ref{fig:SINRvsK}, [Chen2007] has the best performance when very few snapshots are available (e.g., $K<10$).
When $K >  20$, SSC-DL achieves the best performance among all beamforming methods.
The output SINR of SSC-DL appears to be very close to the performance limit depicted in ``Optimal'' when $K>160$, but it is not supposed to approach the limit as $K \to \infty$.
Instead, the curve ``MVDR, no mismatch'' shall approach the limit, as the estimated autocorrelation matrix $\hat{\bf R}_y$ approaches ${\bf R}_y$ as $K\to\infty$.
Therefore, it would eventually exceed the curve ``SSC-DL.''
Although not shown explicitly in Figure \ref{fig:SINRvsK}, our numerical results indicate that this happens only when $K \geq 12,000$.
This shows that the proposed SSC-DL not only is robust to DoA mismatch but also is robust to a limited number of available snapshots.

\subsection{The Projection Ratio of Interference signals}
\label{secsec:ProjectionRatio}
{ In this section, we present the simulation results of projection ratio to validate that the estimated weight vector using the constructed subspace in SSC-DL can more effectively suppress interference than the one estimated by the conventional DL method. The projection ratio for SSC-DL is calculated by the ratio of projections of the steering vector of the interference signal on ${\bf Q}_{s} $ and ${\bf Q}_{n}$:
\begin{align}\label{eq:sscdlpr}
\mbox {Projection Ratio for SSC-DL} =  \frac{\left\| {\bf Q}_{s}{\bf Q}^H_{s}{\bf a}(\theta_{Ij})\right\|}{\left\| {\bf Q}_{n}{\bf Q}^H_{n}{\bf a}(\theta_{Ij})\right\|},
\end{align}
where ${\bf Q}_s$ and ${\bf Q}_n$, as defined in Section \ref{subsec:interf_suppr}, are matrices whose columns form orthonormal bases for the subspaces ${\cal W}$ in SSC-DL and ${\cal W}^\perp$, respectively.} Notably, the projection ratio for SSC-DL indicates the ratio of the closeness of ${\bf a}(\theta_{Ij})$ to ${\bf Q}_s$ over that of ${\bf a}(\theta_{Ij})$ to ${\bf Q}_n$. A smaller projection ratio indicates that ${\bf a}(\theta_{Ij})$ is closer to the matrix ${\bf Q}_n$, which also suggests that ${\bf Q}_s$ of SSC-DL can suppress the interference more effectively.
{Next, to compute the projection ratio for DL, we project the steering vector of interference signal on ${\bf w}_{DL}$ (the weight vector estimated by the DL method) and a matrix ${\bf Q}_{l}$ whose columns form an orthonormal basis for the null space of ${\bf w}_{DL}$. Then the projection ratio for DL is defined by:
\begin{align}\label{eq:dlpr}
\mbox {Projection Ratio for DL} =  \frac{\left\|  \frac{{\bf a}^H(\theta_{Ij}){\bf w}_{DL}}{\left \| {\bf w}_{DL} \right \|^2} {\bf w}_{DL} \right\|}{\left\| {\bf Q}_{l}{\bf Q}^H_{l}{\bf a}(\theta_{Ij})\right\|}.
\end{align}
Similar to the projection ratio for SSC-DL in (\ref{eq:sscdlpr}), the projection ratio for DL indicates the ratio of the closeness of ${\bf a}(\theta_{Ij})$ to ${\bf w}_{DL}$ over that of ${\bf a}(\theta_{Ij})$ to ${\bf Q}_l$. A smaller projection ratio indicates that ${\bf a}(\theta_{Ij})$ is closer to the matrix ${\bf Q}_l$, which also suggests that ${\bf w}_{DL}$ can suppress the interference more effectively.

Figure \ref{fig:prvssnr} and Figure \ref{fig:prvssnapshot} show the results of projection ratios versus SNR(dB) and snapshot number ($K$), respectively. Note that Figure \ref{fig:prvssnapshot} is showed with fixed SNR = 10 dB. The curves of ``Interference $j$-th with the subspace of SSC-DL'' and ``Interference $j$-th with DL'' are computed by (\ref{eq:sscdlpr}) and (\ref{eq:dlpr}), respectively. From these two figures, the curve ``Interference 1 with the subspace of SSC-DL'' yields smaller projection ratios than the curve ``Interference 1 with DL'', confirming that the proposed SSC-DL can estimate the weight vectors in a more effective way to suppress interference 1 than the DL method. Similarly, we note that SSC-DL also provides more suppression for interference 2, again confirming that the proposed SSC-DL outperforms the DL method in terms of suppressing interference signals to achieve better output SINR performance. }
\label{pr}
\begin{figure}[h]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnr_projectionratio}}
\caption{The Projection Ratio of Interference signals versus SNR  }
\label{fig:prvssnr}
\end{figure}
\begin{figure}[h]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnapshot_projectionratio}}
\caption{The Projection Ratio of Interference signals versus snapshot number ($K$).  }
\label{fig:prvssnapshot}
\end{figure}

\subsection{SINR versus subspace dimension}
\label{subsec:SINRvsM}
\begin{figure}[h]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{SINRvsM_snr10}}
\caption{Output SINR versus the subspace dimension $M$.  }
\label{fig:sinrvsm}
\end{figure}

In Figure \ref{fig:sinrvsm}, we study the SINR performance versus the given subspace dimension $M$. As described in Section \ref{subsec:M}, we have learned that $M$ must be large enough to support any possible desired DoA within the range $[\theta_1, \theta_2]$. However, $M$ cannot be too large either: when $M=N$, the SSC-DL degenerates to the conventional DL and has an SINR performance as poor as conventional DL.

\subsection{SINR performance with various DoAs}
\label{subsec:SINRwithVariousDoA}
In this subsection, we display numerical results that show the effectiveness of the proposed SSC-DL method in the environments of various DoAs of the desired signal.

\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnr7}}
\caption{Output SINR versus SNR  when $\theta_d = -45^\circ$. }
\label{fig:snr7}
\end{figure}
\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnapshot7}}
\caption{Output SINR versus snapshot number  when $\theta_d = -45^\circ$. }
\label{fig:snapshot7}
\end{figure}

Figures \ref{fig:snr7} and \ref{fig:snapshot7} show the SINR performance of all methods when the DoA of the desired signal is set to $\theta_d = -45^\circ$, and the DoAs of two interference signals are set to $\theta_{I1} = -21.41^\circ$ and $\theta_{I2} = -11.95^\circ$. The estimated DoA of the desired signal is $\theta_0 = \sin^{-1}\left(\sin\theta_d+\sin 2.5^\circ\right)$ and the values of the boundary angles $\theta_1$ and $\theta_2$ are chosen as $\theta_1 = -36.42^\circ$ and $\theta_2 = -47.16^\circ$, following the formulas
$\theta_1 = \sin^{-1}\left(\sin\theta_0-\sin 4^\circ\right) $ and
$\theta_2 = \sin^{-1}\left(\sin\theta_0+\sin 4^\circ\right).$
The power levels of the interference are chosen the same as described in Section \ref{subsec:simsetup}:
$E[|x_{1}(t)|^{2}]=\sigma_{I_{1}}^{2}=10^4(40\mbox{dB})$ and
$E[|x_{2}(t)|^{2}]=\sigma_{I_{2}}^{2}=10^2(20\mbox{dB})$.

In these plots, we find that the proposed SSC-DL remains its superiority in SINR performance to all other methods.
When the number of snapshots is small ($K<10$), the SINR of SSC-DL is slightly lower than those of [Chen2007] and [Huang2015].
However, when $K>60$, it is observed that the SSC-DL has the best SINR performance among all methods.

\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnr8}}
\caption{Output SINR versus SNR when $\theta_d = -76^\circ$.}
\label{fig:snr8}
\end{figure}

\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnapshot8}}
\caption{Output SINR versus snapshot number when $\theta_d = -76^\circ$.}
\label{fig:snapshot8}
\end{figure}
Figures \ref{fig:snr8} and \ref{fig:snapshot8} display the SINR performance when the desired DoA is $\theta_d = -76^\circ$. This is an angle that a beamformer based on phase arrays would have a poor resolution.
The interference signals have DoAs at $\theta_{I1} = -38.92^\circ$ and $\theta_{I2} = -28.05^\circ$.
The estimated DoA of the desired signal is $\theta_0 = -67.92^\circ$
The range of possible desired signal DoA $[\theta_1, \theta_2]$ where $\theta_1 = -58.97^\circ$ and $\theta_2 = -85.16^\circ$, following the formulas
$\theta_1 = \sin^{-1}\left(\sin\theta_0-\sin 4^\circ\right) $ and
$\theta_2 = \sin^{-1}\left(\sin\theta_0+\sin 4^\circ\right).$
The power levels of the interference are chosen the same as described in Section \ref{subsec:simsetup}:
$E[|x_{1}(t)|^{2}]=\sigma_{I_{1}}^{2}=10^4(40\mbox{dB})$ and
$E[|x_{2}(t)|^{2}]=\sigma_{I_{2}}^{2}=10^2(20\mbox{dB})$.

In the simulation plots, we observe that the proposed method still has the best performance among all methods.
The SINR performance of SSC-DL is very close to that of Optimal as long as a sufficient number of snapshots are available ($K>=200$).
When the number of snapshots is small ($K<10$), the proposed method performs slightly worse than [Huang2015], but when $K>30$, it outperforms all other methods.
{Also, according to simulation results with various settings (not present here due to page limit), we found that generally our method would outperform all other existing methods when $K$ is chosen as $\geq 5 \cdot N$ (where $N$ is the number of antennas); when $K$ is chosen as $\geq 20 \cdot N$, the SINR performance roughly reaches its limit.}

\subsection{SINR performance in a larger signal dimension}

\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnr14}}
\caption{Output SINR versus SNR when $N=14$. }
\label{fig:snr14}
\end{figure}

\begin{figure}[htp]
\centering \centerline{
\includegraphics[width=0.5\textwidth]{sinrvssnapshot14}}
\caption{Output SINR versus Snapshot when $N=14$. }
\label{fig:snapshot14}
\end{figure}

Figures \ref{fig:snr14} and \ref{fig:snapshot14} show the SINR performance when the number of antennas increases to $N = 14$.
The desired DoA is $\theta_d = 0^\circ$. Other simulation parameters are the same as those in Section \ref{subsec:simsetup}.
In this case, the SINR performance of [Chen2007] degrades in the high SNR region since the method \cite{chen2007} requires the parameter $N$ to be upper bounded by a factor in order to have a satisfactory performance.
The proposed SSC-DL still has an excellent performance and outperforms all other methods.

\section{Conclusion} \label{sec:conclusion}
In this paper, we proposed a novel yet simple method, namely subspace-constrained diagonal loading (SSC-DL), to deal with the DoA mismatch issue in a beamforming system.
The proposed method imposes an additional subspace constraint on the conventional diagonal loading (DL) method to provide better interference suppression and has a closed-form solution.
The proposed method maintains a high magnitude response to the desired DoA but does not have to widen the beamwidth like other existing approaches.
The signal-to-interference-plus-noise ratio (SINR) performance, as shown in numerical results, is even approaching the theoretical limit suggested in \cite{capon1969}.
In the future, an analytical study would be helpful in further understanding the superiority of the proposed method.
Further, automatic selection of parameters $M$ and $\gamma$ is challenging yet desirable.
In addition, {making the DoA of the desired signal modeled as a random variable} and extending the proposed SSC-DL method to the case of general steering vector errors is also of interest.

\section*{Statement on Conflict of Interests}
The authors declare that there is no conflict of interests regarding the publication of this paper.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi
\printbibliography

\end{document}