\section{Introduction}\label{introduction}
We are progressing through an exciting period in Quantum Technologies and with the small-scale commercial quantum computers becoming increasingly available \cite{2018arXiv181203050T}, Quantum Information Processing is witnessing spectacular developments.   Before quantum processors become scalable devices capable of error correction and universality \cite{kim2019leveraging}, the current and near-term devices, referred to as the Noisy Intermediate-Scale Quantum (NISQ) \cite{preskill2018quantum} devices are getting explored for solving certain hard problems to achieve significant speedups over the best known classical algorithms \cite{1997quant.ph..1001B}.  Promising results are already reported for solutions in the areas of optimisation, chemistry, machine learning. Apart from only speed up considerations, exploiting quantum mechanical properties of superposition, entanglement and interference for solving problems differently with possible performance improvements are getting explored, among others.  Needless to say, hybrid quantum algorithms which use both classical and quantum resources to solve potentially difficult problems \cite{zhu2019training} are worked out and put into action.

It has been brought out that unsupervised machine learning and the associated optimisation strategies can be elegantly handled by quantum or hybrid quantum algorithms.  In this paper, we consider clustering, an important unsupervised task.  Clustering consists of assigning labels to elements of a dataset based only on how similar they are to each other - like objects will have the same label, unlike objects will have different labels \cite{otterbach2017unsupervised}. In order to represent dissimilarity (or similarity), we need to define a distance measure between  two  data samples.  The distance between every possible  pair  of  data  samples can be captured in a matrix. This matrix can be interpreted as an adjacency matrix of a graph, where each vertex or node represents an element of data set and the the weight of edge between vertices is the corresponding distance \cite{otterbach2017unsupervised}. In clustering, the main assumption is that distant points belong to different clusters; hence maximizing the overall sum of all weights (distances) between nodes with different labels represents a natural clustering algorithm for two-cluster case. The mathematical formulation of this is a well known Maximum-Cut (Max-cut) problem and it can be easily translated to an optimisation objective \cite{otterbach2017unsupervised}. The Max-cut problem  is  an  example  of  the  class of  NP-complete  problems,  which  are  notoriously hard to solve.  The Max-Cut and many other combinatorial problems, like, machine scheduling, computer-aided design, traffic message management \cite{otterbach2017unsupervised} fall under the unifying model of Quadratic Unconstrained Binary Optimisation (QUBO) \cite{otterbach2017unsupervised, lewis2017quadratic}.  One approach to solving Max-cut is to construct a physical system, typically a set of interacting spin  particles (two-state particles) whose lowest energy state encodes the solution to the problem, so that solving the problem is equivalent to finding the ground state of the system.

Two main approaches have been identified to find the ground state of interacting spin systems (quantum optimisation) in NISQs \cite{kim2019leveraging,otterbach2017unsupervised}: Quantum Annealing (QA) and Quantum Approximate Optimisation Algorithms (QAOA) \cite{farhi2014quantum}. QA is a form of analog computation that has been developed theoretically in the early nineties but realized experimentally in a programmable device only in 2011 by D-Wave Systems.  QAOA, invented in 2014 and recently generalised for constrained combinatorial optimisation, requires digital gate-model quantum computing; it can be seen in some parameter range as a “digitised” version of QA \cite{kim2019leveraging}.

In this paper, our starting point is a graph and we examine clustering on this abstraction (the graph itself can be constructed from the data points as cursorily mentioned in the beginning). To start with, we briefly touch upon the 2-cluster max-cut problem in terms of the usual Ising model of interacting spins, but report some additional results/observations related to graph components (independent subgraphs). Then,  we propose a simple way to extend the strategy to address 3-cluster problem on graphs.  The requisite 3-state particles interaction and the associated Hamiltonian are brought out. Apart from providing typical results, remarks on how to go about implementations are also made, including on the hypothetical qutrit computing device.  Extrapolation to more than three cluster case involving qudits is also suggested, culminating in the Quadratic Unconstrianed $D$-ary Optimisation (QUDO).

The paper is organized as follows: In Section \ref{maxcut}, Ising model and the max-cut clustering are presented.  Clustering into $3$ classes is covered in detail in Section \ref{max3cut}. Graph clustering into $d$ classes is brought out in Section \ref{maxkcut}.  In both Section \ref{max3cut} and Section \ref{maxkcut} relevant results are interspersed. Remarks related to implementation are provided in Section  \ref{implementation} including systematic steps for constructing the qudit-based circuit for $d$-ary QAOA to solve max $d$-cut. Conclusions are provided in Section \ref{conclusions}.

\section{Two-Group Clustering Max-Cut Problem}\label{maxcut}
As mentioned in the previous section, one way to solve the two-cluster graph maxcut problem is to have a model of two-state interacting particles and solve for the lowest energy state.  This interaction model (for spins) is the Ising model, originally developed to describe ferromagnetism, but subsequently extended to more problems \cite{1525}.

\subsection{Ising Model}\label{IsingModel}

The Ising model can be formulated on any graph as follows: consider an undirected graph
$G=(V,E)$, where $V=\{v_{1},...,v_{N}\}$  is a set of $N$ sites, and $E$ is a set of edges representing the interactions between these sites. Every site $i$ has a corresponding spin variable $s_{i}$ \cite{kamenetsky2010ising}. These spins are binary-valued, taking values $+1$ for “up” or $-1$ for “down”. Two spins  $s_{i}$ and $s_{j}$  may interact with each other \cite{kamenetsky2010ising}. The energy of such an interaction depends on whether the values of the participating spins are the same or different: it is given by $J_{ij}s_{i}s_{j}$, where $J_{ij}$  is the strength of the interaction \cite{kamenetsky2010ising}.

For each pair of interacting spins $s_{i}$  and $s_{j}$ (i.e., $J_{ij}\neq0$), there exists a corresponding edge $(i,j) \in E$. The state of the model, $\textbf{\emph{s}}$, is an assignment of all $N$ variables $s_{i}$, $1 \leq i \leq N$. The set of all possible configurations is $\xi = \{-1,1\}^N$ \cite{kamenetsky2010ising}. As well as pair-wise interactions, there can also be an external field that affects each site $i$ with energy $h_{i}s_{i}$. Thus, in the general case, the energy of a configuration $\textbf{\emph{s}} \in \xi$  is given by the so-called Edwards-Anderson Hamiltonian \cite{kamenetsky2010ising}:
\begin{equation}
	\label{eq:classicalIsing1}
	H(\textbf{\emph{s}})=\!\!\!\!\sum_{(i,j) \in E}\!\!\!\!J_{ij}s_{i}s_{j} + \sum_{i \in V}h_{i}s_{i}
\end{equation}

When $h_{i} = 0 \; \forall \; i \in V$, the system is said to have no external field (also called zero magnetic field condition), in which case the energy of the configuration $\textbf{\emph{s}}$ becomes:
\begin{equation}
	H(\textbf{\emph{s}})=\!\!\!\!\sum_{(i,j) \in E}\!\!\!\!J_{ij}s_{i}s_{j}
	\label{eq:classicalIsing}
\end{equation}

In this paper, we mostly consider the Hamiltonian with zero external field.  The system prefers lower energy states, i.e., those $\textbf{\emph{s}}$ that minimise $H(\textbf{\emph{s}})$.  An important task is to find configurations that minimise the energy of the system; such a configuration is known as ground state.

In order to obtain the Quantum Mechanical description of the Ising model of Equation \eqref{eq:classicalIsing1} and Equation \eqref{eq:classicalIsing}, one has to replace each $s_{i}$ by the Pauli-$Z$ matrix given by $\sigma_{z}^{i}$. As the state corresponding to the quantum mechanical interaction of one or more particles is given by the tensor products of the corresponding individual states, the $\sigma_{z}^{i} \sigma_{z}^{j}$ terms in $H(\boldsymbol{\sigma})$ denote the tensor product between $\sigma_{z}^{i}$ and $\sigma_{z}^{j}$, where $\boldsymbol{\sigma}$ is the cumulative spin configuration of the complete system. It must be noted that these terms need to be appropriately constructed through tensor products of $\sigma_{z}$ and $I$ matrices, as discussed in subsection \ref{maxcutIsing}, to capture the pair-wise interaction between the $i^{th}$ and $j^{th}$ spins. With these considerations in mind, the Equation \eqref{eq:classicalIsing1} and Equation \eqref{eq:classicalIsing} convert to \cite{atan}:

\begin{equation}
	H(\boldsymbol{\sigma})=\sum_{ij}J_{ij} \sigma_{z}^{i} \sigma_{z}^{j} + \sum_{i}h_{i} \sigma_{z}^{i}
	\label{eq:quantumIsingHamiltonian}
\end{equation}

\begin{equation}
	H(\boldsymbol{\sigma})=\sum_{ij}J_{ij} \sigma_{z}^{i} \sigma_{z}^{j}
	\label{eq:finalQuantumIsingHamiltonian}
\end{equation}

\subsection{Max-cut Problem Using Ising Model}\label{maxcutIsing}
Little more formally, the max-cut problem starts with an undirected graph $G(V,E)$ with a set of vertices $V$ and a set of edges $E$ between the vertices. The weight $w_{ij}$ of an edge between vertices $i$ and $j$ is a positive real number, with $w_{ij} = 0$ if there is no edge between them. A cut is a set of edges that separates the vertices $V$ into two disjoint sets $V_{1}$ and $V_{2}$, such that $V_{1} \subseteq V$ and $V_{2}=V \setminus V_{1}$, and the cost of a cut is defined as the as the sum of all weights of edges connecting vertices in $V_{1}$ with vertices in $V_{2}$. One can connect the cut to the Ising model by identifying the vertices with $s_{i}$ and $w_{ij}$ with $J_{ij}$; $s_{i} = 1$ suggesting that node $i$ belongs to $V_{1}$ and $s_{i} = -1$ corresponds to $V_{2}$ (of course, $V_{1}$ and $V_{2}$ can be interchanged). The cost of the cut can, in fact, be written in terms of the objective function:
\begin{equation}
	C=\sum_{ij}\frac{w_{ij}}{2}(1-s_{i}s_{j})
	\label{eq:maxcutCost}
\end{equation}
The max-cut problem aims at partitioning the nodes in such a way that the cost of the resulting cut is maximised. On comparing Equation \eqref{eq:classicalIsing} and Equation \eqref{eq:maxcutCost}, it must be noted that the sign of the $s_{i}s_{j}$ term changes. Thus, the max-cut then boils down to finding the lowest energy state, instead of the maximum, of Ising model with $J_{ij} = w_{ij}$:

\begin{equation}
	\max C = \min H(\textbf{\emph{s}})
\end{equation}

\noindent by noting that scaling the cost function by a constant multiplicative factor does not change the optimal solutions. As remarked earlier, the Max-Cut problem is equivalent to QUBO formulation where the two-state variable is $0$ or $1$ binary valued instead of $1$ and $-1$; if the QUBO variables are denoted by $x_{i}$, the two formulations are related by $s_{i}=2x_{i}-1$ \cite{lewis2017quadratic,ushijima2017graph,gowda}.

The classical Ising max-cut can be extended to the quantum framework by replacing $s_{i}$ by $\sigma_{z}^{i}$, as discussed in Section \ref{IsingModel}. This gives us the final quantum Ising Hamiltonian, given by Equation \eqref{eq:IsingHamiltonian}, which needs to be minimised to get the optimum cut \cite{ushijima2017graph}. In carrying out the optimisation based on Equation \eqref{eq:IsingHamiltonian}, if two adjacent nodes fall into the same cluster, then $w_{ij}$ is added to the cost function, else it is subtracted.
\begin{equation}
	\label{eq:IsingHamiltonian}
	H(\boldsymbol{\sigma})=\sum_{ij}w_{ij} \sigma_{z}^{i} \sigma_{z}^{j}
\end{equation}

As discussed in the previous subsection, the $\sigma_{z}^{i} \sigma_{z}^{j}$ terms represent the interaction of the nodes through the use of tensor products, as shown in the following example. Let us consider a graph having $5$ vertices. If there is an edge between the nodes $v_{1}$ and $v_{3}$ having weight $w_{13}$, then the interaction between the respective vertices is given by:
\begin{equation}
	H_{13}=w_{13}(I \otimes \sigma_{z} \otimes I \otimes \sigma_{z} \otimes I )
\end{equation}
The Hamiltonian thus formed is always diagonal, as there are no $\sigma_{x}$ terms involved, owing to the tensor product between diagonal matrices. It is to be noted that Ising models can be of transverse field type, where, $\sigma_x$ terms corresponding to Pauli-X matrix also are present, in which case the corresponding Hamiltonian would be non diagonal.

\begin{figure}
	\centering
	\includegraphics[scale=0.59]{binaryOutput}
	\caption{Result of max-cut clustering on weighted graphs. The colour of each node denotes the cluster it belongs to.}
	\label{fig:binary}
\end{figure}

The cost Hamiltonian $H(\boldsymbol{\sigma})$ of the graph $G$, being diagonal, has orthogonal eigenvectors that form a complete standard basis. The system settles to the state having the lowest energy and the eigenvector corresponding to it is given by a unit vector along a standard basis state. The least-energy eigenvector is a $2^{N} \times 1$ vector, whose ket representation gives us an $N$-length bit-string. The nodes of the graph, $G$, are labelled as $0$ or $1$, according to the digits in the bit-string, with the most significant bit representing the label of the first vertex. This gives us a binary-clustered graph. The result of partitioning a graph with the given algorithm has been shown in Figure \ref{fig:binary}. The minimum energy state corresponding to the partition is given by $\ket{010101100101}$.

\section{Solving Max 3-cut using Ising Model}\label{max3cut}
In the direction of arriving at the quantum-assisted solution for Max $3$-cut problem, few observations are put in place. Consider the nodes $2$ and $3$ of a graph with $6$ nodes, shown in Figure \ref{fig:2nodes} for binary clustering:

\begin{figure}[h!]
	\centering
	\includegraphics[scale=0.5]{2nodes}
	\caption{Two nodes connected by an edge}
	\label{fig:2nodes}
\end{figure}

The cost for the connection, in the Ising model, is given by:
\begin{equation}
	H_{23}=5 \times (I \otimes I \otimes \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \otimes \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} \otimes I \otimes I)
\end{equation}
which can be simplified as:
\begin{equation}
	H_{23}=5 \times (I \otimes I \otimes \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \otimes I \otimes I)
\end{equation}
The $4 \times 4$ matrix, in the equation above, incorporates all the possible cluster combinations of the nodes 2 and 3 (along the diagonal). If the two cluster are named $0$ and $1$, then their combinations and energies are shown in Figure \ref{fig:binaryHamiltonian}.

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.5]{binaryHamiltonian.png}
	\caption{Node-node interactions and corresponding interaction energies}
	\label{fig:binaryHamiltonian}
\end{figure}
The energy is $1$ when both the nodes are classified into the same cluster, $\ket{00}$ or $\ket{11}$, and is $-1$ otherwise. The system will settle for the lower energy state, $-1$, thus providing the optimum clustering.

The idea can be extended to clustering the nodes of a graph into $3$ classes, with class labels $0$, $1$ and $2$. The interaction-matrix between adjacent nodes of the graph should ideally look like (coupling strength has been considered to be unity) the matrix in Figure \ref{fig:ternaryHamiltonian}.

\begin{figure}[h]
	\centering
	\includegraphics[scale=0.65]{ternaryHamiltonian.png}
	\caption{Node-node interaction energies for the ternary clustering case}
	\label{fig:ternaryHamiltonian}
\end{figure}

In the binary case, the nodes were represented by $2 \times 2$ Pauli-$Z$ matrices, the tensor product of which gave the required Hamiltonian matrix. Thus, for the 3-class problem, we need $3 \times 3$ matrices for each node, whose tensor product with another such matrix can give the requisite $9 \times 9$ Hamiltonian.

One can think of using a matrix with the cube roots of unity placed along the diagonal of the $3 \times 3$ matrix:
\begin{equation}
	\Omega_{3}=\begin{bmatrix} 1 & 0 & 0 \\ 0 & e^{\nicefrac{2 \pi i}{3}} & 0 \\ 0 & 0 & e^{\nicefrac{4 \pi i}{3}} \end{bmatrix}
\end{equation}

In an interaction between two nodes, the first node is represented by $\Omega_{3}$ and the second node is represented by the complex conjugate transpose of $\Omega_{3}$ , $\Omega_{3}^{\dagger}$ . The 3 in subscript signifies that the cube roots of unity are used to form the matrix. Thus, the energy of interaction of nodes 2 and 3 is:
\begin{equation}
	H_{23}=\dots \otimes \begin{bmatrix} 1 & 0 & 0 \\ 0 & e^{\nicefrac{2 \pi i}{3}} & 0 \\ 0 & 0 & e^{\nicefrac{4 \pi i}{3}} \end{bmatrix} \otimes \begin{bmatrix} 1 & 0 & 0 \\ 0 & e^{\nicefrac{-2 \pi i}{3}} & 0 \\ 0 & 0 & e^{\nicefrac{-4 \pi i}{3}} \end{bmatrix} \otimes \dots
\end{equation}
\begin{equation}
	\Rightarrow H_{23}=\dots \otimes \begin{bmatrix} 1 & 0 & \dots & 0 & 0 \\ 0 & e^{\nicefrac{-2 \pi i}{3}} & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & e^{\nicefrac{2 \pi i}{3}} & 0 \\0 & 0 & \dots & 0 & 1 \\ \end{bmatrix} \otimes \dots
\end{equation}

But, the Hamiltonian, being an observable, must be Hermitian. This can be taken care of by modelling the interaction between two nodes as the tensor product between $\Omega_{3}$ and $\Omega_{3}^{\dagger}$ and taking only the real part of the elements of the resultant matrix. The final form of the interaction is given as:
\begin{equation}
	H_{23}=5 \times I \otimes I \otimes \frac{1}{2} (\Omega_{3}^{2} \otimes \Omega_{3}^{3 \dagger} + \Omega_{3}^{2 \dagger} \otimes \Omega_{3}^{3}) \otimes I \otimes I
\end{equation}

It is interesting to note that this is completely analogous to the $d=2$ case, where the Pauli-$Z$ matrix is Hermitian and thus, $\sigma_{z} \otimes \sigma_{z}^{\dagger} = \sigma_{z}^{\dagger} \otimes \sigma_{z}$.

The term $\bar{H}_{23}= \frac{1}{2} (\Omega_{3}^{2} \otimes \Omega_{3}^{3 \dagger} + \Omega_{3}^{2 \dagger} \otimes \Omega_{3}^{3})$ evaluates to:
\begin{equation}
	\bar{H}_{23}=\begin{bmatrix} 1 & 0 & \dots & 0 & 0 \\ 0 & -0.5 & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & -0.5 & 0 \\0 & 0 & \dots & 0 & 1 \\ \end{bmatrix}
\end{equation}

If adjacent nodes are placed in dissimilar clusters, the interaction energy is $-0.5$, and $1$ for similar clusters. This is in contrast to the desirable energy values of $-1$ and $1$ for dissimilar and similar clusters, respectively. But this does not have any effect on the clustering result since the energy for adjacent vertices in different clusters is still lower than that for similar clusters, and the former will be energetically favoured. To provide an easy visualisation, the cluster classes can be modelled along the cube roots of unity as shown in Figure \ref{fig:ternaryVector}, since the clustering energy can be considered in terms of angular separation, as discussed in the next section.

\begin{figure}
	\centering
	\includegraphics[scale=0.65]{ternaryVector}
	\caption{The cluster classes are modelled along the vectors: $0$ along red, $1$ along blue, $2$ along green.}
	\label{fig:ternaryVector}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[scale=0.5]{ternaryOutput}
	\caption{Results of max $3$-cut clustering where the three clusters are denoted by red, blue and green colours.}
	\label{fig:ternary}
\end{figure}

The typical result of max 3-cut partitioning on a graph can be seen in Figure \ref{fig:ternary}. The final Hamiltonian for max-cut is given by:

\begin{equation}
	\label{eq:max3cutHamiltonian}
	H=\sum_{ij}\frac{w_{ij}}{2}(\Omega_{3}^{i} \Omega_{3}^{j \dagger}+\Omega_{3}^{i \dagger} \Omega_{3}^{j})
\end{equation}

Extensive numerical evaluation studies with different graphs corroborates the applicability of the proposal.

\section{Max $d$-Cut}\label{maxkcut}

Having looked at max-cut and max $3$-cut, where the nodes of a graph are partitioned into two and three classes respectively, in this section, we propose a way of generalising the max-cut problem for $d$ classes. Such problems have traditionally been known as max $d$-cut \cite{Alan}. Before moving onto the problem, it is necessary to introduce some additional mathematical machinery to facilitate the solution.

In section \ref{max3cut}, the matrix $\Omega_{3}$ can be recognised as the clock matrix for three dimensions. This is not surprising as the clock and shift matrices, $U_{d}$ and $V_{d}$, have been used to generalise the Pauli-$Z$ and $X$ matrices \cite{prakash2018normal}, respectively, for $d$ dimensional \emph{qudits} \cite{Frydryszak_2017}. These matrices, which are zero-trace and symmetrical, but not Hermitian, can be utilised to solve $d$-ary optimisation problems where the variables can take $d$ number of discrete values, similar to a $d$-level system. In this section we extend the idea of Quadratic Unconstrained Binary Optimisation, $QU\!BO$, to such $d$-ary problems, and dub the technique as Quadratic Unconstrained $D$-ary Optimisation, $QU\!DO$. To construct the Hamiltonian for the max $d$-cut problem, the $d$-dimensional Clock matrices (specified in Equation \eqref{eq:clock}  \cite{anwar2012qutrit}) are used to appropriately replace the Pauli-$Z$ matrix in Equation \eqref{eq:IsingHamiltonian}, as was done for Max 3-Cut in Equation \eqref{eq:max3cutHamiltonian}.

\begin{equation}
	V_{d}=\sum_{j=0}^{d-1}\ket{j}\bra{(j+1) \; mod \; d}
	\label{eq:shift}
\end{equation}
\begin{equation}
	U_{d}=\sum_{j=0}^{d-1}\omega^{j}\ket{j}\bra{j}
	\label{eq:clock}
\end{equation}
where  $\omega=e^{\nicefrac{2\pi i}{d}}$ is the $d^{th}$ root of unity \cite{anwar2012qutrit}.

The resultant Ising Hamiltonians for max $d$-cut is specified as:
\begin{equation}
	H=\sum_{ij}\frac{w_{ij}}{2}(U_{d}^{i}U_{d}^{j \dagger}+U_{d}^{i \dagger}U_{d}^{j})
	\label{eq:maxkcut}
\end{equation}

Further, it is possible to view the proposed method as a quantum mechanical extension of the vector Potts model \cite{Grimmett2006,wu1982potts}, where the spin-states of a $d$ level system are equally spaced on a unit circle, with the interaction energy of adjacent spins depending on the cosine of the relative angles between the states \cite{wu1982potts}. When the Ising Model contains \textit{external, longitudinal magnetic field} tems, Equantions \eqref{eq:IsingHamiltonian} and \eqref{eq:maxkcut} modify to Equations \eqref{eq:ising} and \eqref{eq:isingd}, respectively, where $h_i$ can be considered as an additional weight of the $i^{th}$ node of the graph under consideration. This is carried out by leveraging Equation \eqref{eq:1sigma}.

\begin{equation}
	\sigma_{z} = \frac{1}{2}(\sigma_{z} + \sigma_{z}^\dagger)
	\label{eq:1sigma}
\end{equation}

\begin{equation}
	H(\boldsymbol{\sigma})=\sum_{ij}w_{ij} \sigma_{z}^{i} \sigma_{z}^{j} + \sum_i h_i \sigma_{z}^i
	\label{eq:ising}
\end{equation}

\begin{equation}
	H=\sum_{ij}\frac{w_{ij}}{2}(U_{d}^{i}U_{d}^{j \dagger}+U_{d}^{i \dagger}U_{d}^{j}) + \sum_i \frac{h_i}{2}(U_d^i + U_d^{i \dagger})
	\label{eq:isingd}
\end{equation}

\begin{figure}
	\centering
	\includegraphics[scale=0.5]{quaternaryOutput2}
	\caption{Results of max $4$-cut clustering on an unweighted graph where the four clusters are denoted by red, blue, green and cyan colours.}
	\label{fig:quaternary2}
\end{figure}

\begin{figure*}
	\centering
	\includegraphics[scale=0.63]{quaternaryHamiltonian}
	\caption{Node-node interaction for the quaternary clustering case}
	\label{fig:quaternaryHamiltonian}
\end{figure*}

A typical result of quaternary clustering for a simple graph is shown in Figure \ref{fig:quaternary2}. Again, extensive verification of the results of the formulation using different graphs and weights has been carried out.

For the $d=4$ case, the interaction Hamiltonian matrix for two adjacent nodes has been shown in Figure \ref{fig:quaternaryHamiltonian}. It should be noted that the $\ket{aa}$ elements are $1$, while the $\ket{ab}$ elements are $0$ or $-1$. This happens because there are two possible angles between the $4^{th}$ roots of unity, i.e., $\nicefrac{\pi}{2}$ and $\pi$. If the angle between the classes is $\nicefrac{\pi}{2}$, then the interaction term is 0, and it is -1 if the classes are $\pi$ radians apart. This means that having a larger angular difference between the classes is more favorable. The effect of clustering still remains the same, however. The system settles for a state that ensures the highest angular difference between the classes. This was not apparent for the $d=2$ or $d=3$ cases because there was only one possible angle between the classes.

For a $d$-cluster problem, if the classes are numbered from $0$ to $(d-1)$, then the interaction energy term between nodes of classes $a$ and $b$ is given by:
\begin{equation}
	z_{a} \cdot z_{b}=\frac{1}{2}(\bar{z_{a}}z_{b}+z_{a}\bar{z_{b}})
	\label{dotproduct}
\end{equation}

where $z_{a}=e^{\nicefrac{2\pi ai}{d}}$ and $\bar{z_{a}}$ is the complex-conjugate of $z_{a}$.

\begin{figure}
	\centering
	\includegraphics[scale=0.65]{quaternaryVector}
	\caption{The cluster classes along the vectors: $0$ along red, $1$ along blue, $2$ along green, and $3$ along cyan.}
	\label{fig:quaternaryVector}
\end{figure}

Figure \ref{fig:quaternaryVector} shows the vector representation of the four classes along the fourth roots of unity.

Additionally, for $d=4$, an interesting behaviour was observed.  For many graphs, the result was bipartite, i.e., the resultant partition had only two classes. The partitioning for such a graph has been shown in Figure \ref{fig:quaternary1}. Other graphs had solution states for $2$, $3$ and $4$ classes, all having the minimum energy eigenvalue.
This might have been due to the fact that adjacent classes have $0$ interaction energy. Further investigations are necessary for more than $4$ classes and to arrive at the nature of the behaviour for general $d$ classes.

\begin{figure}
	\centering
	\includegraphics[scale=0.5]{quaternaryOutput1}
	\caption{Results of max $4$-cut clustering where the nodes get placed into two clusters only, instead of four.}
	\label{fig:quaternary1}
\end{figure}

\section{Proposition for Future Implentation of QUDO Algorithms} \label{implementation}
With the NISQ devices already being available, it is really essential to consider the implementation of the quantum algorithms on the hardware. Few pointers are put forth in this section. The binary max-cut problem can be solved using quantum annealing such that the system settles to the final Hamiltonian given by the Ising function. The D-Wave annealer Hamiltonian may be represented as \cite{dwave}:

\begin{equation}
	H_{ising}=-\frac{A(s)}{2}\left(\sum \sigma_{x}^{i}\right)+\frac{B(s)}{2}\left(\sum_{i} h_{i}\sigma_{z}^{i}+\sum_{i>j}J_{ij}\sigma_{z}^{i}\sigma_{z}^{j}\right)
	\label{eq:dwave}
\end{equation}

The cost function can easily be mapped to the D-Wave Ising Hamiltonian and the minimum energy states can be obtained by appropriate sampling.

Alternatively, the QAOA approach can also be utilised to find the solution of binary max-cut problem by evolving the appropriately prepared wavefunction using unitary operators towards the Hamiltonian obtained by casting the problem into Ising model (or the equivalent QUBO) \cite{hadfield2018quantum,zhou2018quantum}. Thus, QA and QAOA have formulation of the Ising model in common but they require different hardware. Ruslan Shaydulin, et al, have compared the two approaches for a different, but similar, problems \cite{shaydulin2018community}. In QAOA hybrid quantum and classical processing is utilised and the quantum computation is in terms of the gate-circuit model \cite{hadfield2018quantum,cervera2018exact}.

For max $d$-cut (or QUDO) problems, these approaches cannot be used directly. The existing
quantum annealers are inherently binary in nature, due to the presence of $2 \times 2$ Pauli matrices.
Ushijima-Mwesigwa, et al \cite{ushijima2017graph}, have proposed a concept of super nodes for Graph Partitioning into $d$
classes. The problem formulation is quite similar to that of max-cut, and the same
approach can be used for the latter. The drawback is that for a graph with $N$ vertices, $dN$ qubits
are required to model the problem onto an annealer and the matrix blows up by a factor of $d^{2}$.

If useful, future annealers can be designed in such a way that they utilise Equation \eqref{eq:maxkcut} as the
final Hamiltonian, with the addition of cross terms containing $U_{d}$ and $V_{d}$. Specifically, one can
think of annealers based on \emph{qutrits} \cite{Burlakov2003} to address max $3$-cut based on Equation \eqref{eq:max3cutHamiltonian}.

\subsection{Qudit-based Circuit of QAOA for Max $d$-Cut}
In this section, we elaborate on the steps to construct the gate-based quantum circuit for solving the problem of Max $d$-Cut using QAOA \cite{farhi2014quantum}, which utilizes a hybrid approach by leveraging the variational principle, for the $d$-ary case. In QAOA a $p$-layer ansatz $\ket{\psi(\boldsymbol{\gamma}, \boldsymbol{\beta})}$, with the cost and mixing Hamiltonians $H_c$ and $H_x$ are used, where:
\begin{equation}
	\ket{	\psi(\boldsymbol{\gamma}, \boldsymbol{\beta}) }= e^{-i \beta_p H_x}e^{-i \gamma_p H_c} \dots e^{-i \beta_1 H_x}e^{-i \gamma_1 H_c} \ket{+}^{\otimes n}
\end{equation}

The parameters $\boldsymbol\gamma$ and $\boldsymbol\beta$ are found using classical optimizers, subject to the minimization of $\braket{\psi(\boldsymbol{\gamma}, \boldsymbol{\beta})| H_c |\psi(\boldsymbol{\gamma}, \boldsymbol{\beta})}$. The modifications required to use QAOA on qudit-based systems to solve the Max $d$-Cut problem will be described in the following subsections.

\subsubsection{Preparation of the state $\ket{+}_d^{\otimes n}$}

It is widely known that $\ket{+}^{\otimes n}$ is used to refer to the zero-phase, equal-superposition state of the system, and is achievable by the application of the $H$ gate on all the qubits which are initially in the $\ket{0}$ state, i.e.:
\begin{equation}
	\ket{+}^{\otimes n} = H^{\otimes n} \ket{0}^{\otimes n}
\end{equation}

To obtain a similar superposition for qudits, we propose the use of the generalized Walsh-Hadamard matrix $W_d$ \cite{sylvester}, in a $d$-dimensional Hilbert space, where $\omega = e^{\sfrac{2\pi i}{d}}$. To distinguish qudits from the traditionally used qubits, we represent the qudit states as $\ket{x}_d$. Figure \ref{fig:qudit_superposition} shows the circuit to obtain the state $\ket{+}_d$ from $\ket{0}_d$.

\begin{equation}
	W_d = \frac{1}{\sqrt{d}}
	\begin{bmatrix}
		1 & 1 & 1 &\dots& 1 \\
		1 & \omega^{d-1} & \omega^{2(d-1)} &\dots& \omega^{(d-1)^2}\\
		1 & \omega^{d-2} & \omega^{2(d-2)} &\dots& \omega^{(d-1)(d-2)}\\
		\vdots & \vdots & \vdots &\ddots& \vdots \\
		1 & \omega & \omega^2 & \dots & \omega^{d-1} \\
	\end{bmatrix}
\end{equation}

\begin{equation}
	\ket{+}_d^{\otimes n} = W_d^{\otimes n} \ket{0}_d^{\otimes n}
\end{equation}

\begin{figure}[H]
	\centering
	\includegraphics[scale=0.65]{qudit_superposition}
	\caption{Circuit to place qudits in equal, zero-phase superposition from $\ket{0}_d$.}
	\label{fig:qudit_superposition}
\end{figure}

\subsubsection{Construction of $e^{-i \gamma_m H_c}$}

For the given max $d$-cut unitary to be implemented using qudit-circuits, it is imperative to define the following new unitary gates:

A two-qudit controlled-$V_d$ gate that applies the $V_d^k$ gate on the target-qudit when the control-qudit is in the state $|k\rangle_d$. This gate is analogous to the traditional CX (Controlled-NOT) gate in the qubit case. The matrix-form of the gate, with the $i^{th}$ qudit as control, and $j^{th}$ qudit as the target, is given by $CV_d^{(ij)} = diag(I,V_d,V_d^2,\dots,V_d^{d-1})$, where $V_d = \sum_{l=0}^{d-1} |l \rangle \langle (l+1) \; mod \; d|$ and $V_d^k = \sum_{l=0}^{d-1} |l \rangle \langle (l+k) \; mod \; d|$, which extends from Equation \eqref{eq:shift}.

The other gate $\mathcal{U}_d$ takes $d$ parameters and is the qudit counterpart of the $R_z$ gate; $R_z(\theta)$ is a single-qubit rotation through angle $\theta$ (radians) around the $z$-axis. If the parameters are given by $\Phi = \{\phi_0, \phi_1, \dots, \phi_{d-1} \}$, then $\mathcal{U}_d = diag(e^{i\phi_0}, e^{i\phi_1}, \dots, e^{i\phi_{d-1}})$. It may immediately be observed that $\mathcal{U}_d = U_d$, when $\Phi = \{0, \frac{2\pi}{d}, \frac{3\pi}{d}, \dots, \frac{2(d-1)\pi}{d} \}$.

In the $d$-ary case for max-cut, discussed in Section \ref{maxkcut}, each node in the graph is represented by a qudit in the quantum circuit. In the venture to solve max $d$-cut using QAOA, the Ising Hamiltonian derived in Equation \eqref{maxkcut} functions as $H_c$, from which the evolution operator  $e^{i\gamma_m H_c}$ is constructed as shown in Equation \eqref{eq:qaoa_oracle}. The corresponding circuit for a single edge between the nodes $0$ and $1$ of a graph has also been portrayed in Figure \ref{fig:qudit_oracle} with parameter $\gamma_m$ for the $m^{th}$ layer of QAOA. The same argument can be extended to form the cost evolution operator for every edge in the graph and each layer.
\begin{gather}
	e^{-i \gamma_m H_c} =
	e^{-i \gamma_m \sum_{ij} \frac{w_{ij}}{2}(U_dU_d^\dagger+U_d^\dagger U_d)} = \nonumber \\
	\prod_{ij} CV_d^{(ij)} \cdot I \otimes \mathcal{U}_{d\Phi} \cdot CV_d^{\dagger (ij)}
	\label{eq:qaoa_oracle}
\end{gather}
where \\
$\Phi = \{-\gamma_m w_{ij}, -\gamma_m w_{ij} cos(\sfrac{2\pi}{d}), \dots, -\gamma_m w_{ij} cos(\sfrac{2(d-1)\pi}{d})\} $.

\begin{figure*}
	\centering
	\includegraphics[scale=0.6]{qudit_oracle}
	\caption{Circuit to apply the $e^{-i \gamma_m H_c}$ operator  on qudits $0$ and $1$ for an edge between the corresponding nodes.}
	\label{fig:qudit_oracle}
\end{figure*}

\subsubsection{Construction of Mixing Operator, $e^{-i \beta_m H_x}$}

QAOA, being a trotterized, adiabatic annealing procedure, relies on the use of a mixing or driver Hamiltonian $H_x$. Traditionally, the single-qubit version of Grover's Diffusion Operator has been used as the mixing operator. This is evident from the fact that $e^{i\sfrac{\beta X}{2}}= R_x(-\beta) = H((1-e^{i\beta})|0\rangle \langle0| - I)H$, upto a global phase factor \cite{akshay};  where $R_x (\theta)$ signifies a rotation of the state about the $x$-axis by an angle of $\theta$ radians. With various extensions of the QAOA algorithm, better mixing Hamiltonians -- which are more suited to the problem at hand -- have been proposed \cite{quantum_alternate}. However, in this implementation, we will consider the single-qudit version of diffusion operator:
\begin{equation}
	e^{-i\beta_m H_x} = \prod_j W_d \cdot \mathcal{U}_{d(-2\beta_m, 0, \dots, 0)} \cdot W_d^\dagger
\end{equation}

\begin{figure}[H]
	\centering
	\includegraphics[scale=0.45]{qudit_mixing}
	\caption{Circuit to implement the $e^{-i \beta_m H_x}$ operator}
	\label{fig:qudit_mixing}
\end{figure}

The circuit in Figure \ref{fig:qudit_mixing} shows the procedure to apply the mixing operator in the $m^{th}$ layer to two of the qudits with parameter $\beta_m$. Again, the same circuit with corresponding parameters can be used for ach of the $p$ layers in QAOA. The absence of qudit-based quantum processors and simulators has deterred us from practically trying out the circuits. However, as mentioned in earlier sections, the numerical evaluation of the results based on the suggested formulation did bring out the usefulness of the approach.

This concludes the discussion related to both annealer and gate-model based implementation suggestions for QUDO formalism.

\section{Conclusions}\label{conclusions}

The paper, by appropriately combining mathematical aspects and supporting numerical results proposed quantum-assisted graph clustering for three or more clusters. Even though the presentation is biased towards algorithmic aspects, pointers are provided for possible architectures for implementation. In the process, qudit based circuit for max $d$-cut pivoted on Quantum Approximate Optimization Algorithm is worked out. While firm footing has been established for max $3$-cut case, for more than $3$ clusters one can carry out additional research to refine the clustering results.

\section*{Acknowledgment}

The authors sincerely thank Mr. Mahesh Rangarajan, Dr. Arpan Pal and Dr. Balamuralidhar P of TCS R\&I for their support and encouragement.

\section{Introduction}
This document is a model and instructions for \LaTeX.
Please observe the conference page limits.

\section{Ease of Use}

\subsection{Maintaining the Integrity of the Specifications}

The IEEEtran class file is used to format your paper and style the text. All margins,
column widths, line spaces, and text fonts are prescribed; please do not
alter them. You may note peculiarities. For example, the head margin
measures proportionately more than is customary. This measurement
and others are deliberate, using specifications that anticipate your paper
as one part of the entire proceedings, and not as an independent document.
Please do not revise any of the current designations.

\section{Prepare Your Paper Before Styling}
Before you begin to format your paper, first write and save the content as a
separate text file. Complete all content and organizational editing before
formatting. Please note sections \ref{AA}--\ref{SCM} below for more information on
proofreading, spelling and grammar.

Keep your text and graphic files separate until after the text has been
formatted and styled. Do not number text heads---{\LaTeX} will do that
for you.

\subsection{Abbreviations and Acronyms}\label{AA}
Define abbreviations and acronyms the first time they are used in the text,
even after they have been defined in the abstract. Abbreviations such as
IEEE, SI, MKS, CGS, ac, dc, and rms do not have to be defined. Do not use
abbreviations in the title or heads unless they are unavoidable.

\subsection{Units}
\begin{itemize}
\item Use either SI (MKS) or CGS as primary units. (SI units are encouraged.) English units may be used as secondary units (in parentheses). An exception would be the use of English units as identifiers in trade, such as ``3.5-inch disk drive''.
\item Avoid combining SI and CGS units, such as current in amperes and magnetic field in oersteds. This often leads to confusion because equations do not balance dimensionally. If you must use mixed units, clearly state the units for each quantity that you use in an equation.
\item Do not mix complete spellings and abbreviations of units: ``Wb/m\textsuperscript{2}'' or ``webers per square meter'', not ``webers/m\textsuperscript{2}''. Spell out units when they appear in text: ``. . . a few henries'', not ``. . . a few H''.
\item Use a zero before decimal points: ``0.25'', not ``.25''. Use ``cm\textsuperscript{3}'', not ``cc''.)
\end{itemize}

\subsection{Equations}
Number equations consecutively. To make your
equations more compact, you may use the solidus (~/~), the exp function, or
appropriate exponents. Italicize Roman symbols for quantities and variables,
but not Greek symbols. Use a long dash rather than a hyphen for a minus
sign. Punctuate equations with commas or periods when they are part of a
sentence, as in:
\begin{equation}
a+b=\gamma\label{eq}
\end{equation}

Be sure that the
symbols in your equation have been defined before or immediately following
the equation. Use ``\eqref{eq}'', not ``Eq.~\eqref{eq}'' or ``equation \eqref{eq}'', except at
the beginning of a sentence: ``Equation \eqref{eq} is . . .''

\subsection{\LaTeX-Specific Advice}

Please use ``soft'' (e.g., \verb|\eqref{Eq}|) cross references instead
of ``hard'' references (e.g., \verb|(1)|). That will make it possible
to combine sections, add equations, or change the order of figures or
citations without having to go through the file line by line.

Please don't use the \verb|{eqnarray}| equation environment. Use
\verb|{align}| or \verb|{IEEEeqnarray}| instead. The \verb|{eqnarray}|
environment leaves unsightly spaces around relation symbols.

Please note that the \verb|{subequations}| environment in {\LaTeX}
will increment the main equation counter even when there are no
equation numbers displayed. If you forget that, you might write an
article in which the equation numbers skip from (17) to (20), causing
the copy editors to wonder if you've discovered a new method of
counting.

{\BibTeX} does not work by magic. It doesn't get the bibliographic
data from thin air but from .bib files. If you use {\BibTeX} to produce a
bibliography you must send the .bib files.

{\LaTeX} can't read your mind. If you assign the same label to a
subsubsection and a table, you might find that Table I has been cross
referenced as Table IV-B3.

{\LaTeX} does not have precognitive abilities. If you put a
\verb|\label| command before the command that updates the counter it's
supposed to be using, the label will pick up the last counter to be
cross referenced instead. In particular, a \verb|\label| command
should not go before the caption of a figure or a table.

Do not use \verb|\nonumber| inside the \verb|{array}| environment. It
will not stop equation numbers inside \verb|{array}| (there won't be
any anyway) and it might stop a wanted equation number in the
surrounding equation.

\subsection{Some Common Mistakes}\label{SCM}
\begin{itemize}
\item The word ``data'' is plural, not singular.
\item The subscript for the permeability of vacuum $\mu_{0}$, and other common scientific constants, is zero with subscript formatting, not a lowercase letter ``o''.
\item In American English, commas, semicolons, periods, question and exclamation marks are located within quotation marks only when a complete thought or name is cited, such as a title or full quotation. When quotation marks are used, instead of a bold or italic typeface, to highlight a word or phrase, punctuation should appear outside of the quotation marks. A parenthetical phrase or statement at the end of a sentence is punctuated outside of the closing parenthesis (like this). (A parenthetical sentence is punctuated within the parentheses.)
\item A graph within a graph is an ``inset'', not an ``insert''. The word alternatively is preferred to the word ``alternately'' (unless you really mean something that alternates).
\item Do not use the word ``essentially'' to mean ``approximately'' or ``effectively''.
\item In your paper title, if the words ``that uses'' can accurately replace the word ``using'', capitalize the ``u''; if not, keep using lower-cased.
\item Be aware of the different meanings of the homophones ``affect'' and ``effect'', ``complement'' and ``compliment'', ``discreet'' and ``discrete'', ``principal'' and ``principle''.
\item Do not confuse ``imply'' and ``infer''.
\item The prefix ``non'' is not a word; it should be joined to the word it modifies, usually without a hyphen.
\item There is no period after the ``et'' in the Latin abbreviation ``et al.''.
\item The abbreviation ``i.e.'' means ``that is'', and the abbreviation ``e.g.'' means ``for example''.
\end{itemize}
An excellent style manual for science writers is \cite{b7}.

\subsection{Authors and Affiliations}
\textbf{The class file is designed for, but not limited to, six authors.} A
minimum of one author is required for all conference articles. Author names
should be listed starting from left to right and then moving down to the
next line. This is the author sequence that will be used in future citations
and by indexing services. Names should not be listed in columns nor group by
affiliation. Please keep your affiliations as succinct as possible (for
example, do not differentiate among departments of the same organization).

\subsection{Identify the Headings}
Headings, or heads, are organizational devices that guide the reader through
your paper. There are two types: component heads and text heads.

Component heads identify the different components of your paper and are not
topically subordinate to each other. Examples include Acknowledgments and
References and, for these, the correct style to use is ``Heading 5''. Use
``figure caption'' for your Figure captions, and ``table head'' for your
table title. Run-in heads, such as ``Abstract'', will require you to apply a
style (in this case, italic) in addition to the style provided by the drop
down menu to differentiate the head from the text.

Text heads organize the topics on a relational, hierarchical basis. For
example, the paper title is the primary text head because all subsequent
material relates and elaborates on this one topic. If there are two or more
sub-topics, the next level head (uppercase Roman numerals) should be used
and, conversely, if there are not at least two sub-topics, then no subheads
should be introduced.

\subsection{Figures and Tables}
\paragraph{Positioning Figures and Tables} Place figures and tables at the top and
bottom of columns. Avoid placing them in the middle of columns. Large
figures and tables may span across both columns. Figure captions should be
below the figures; table heads should appear above the tables. Insert
figures and tables after they are cited in the text. Use the abbreviation
``Fig.~\ref{fig}'', even at the beginning of a sentence.

\begin{table}[htbp]
\caption{Table Type Styles}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Table}&\multicolumn{3}{|c|}{\textbf{Table Column Head}} \\
\cline{2-4}
\textbf{Head} & \textbf{\textit{Table column subhead}}& \textbf{\textit{Subhead}}& \textbf{\textit{Subhead}} \\
\hline
copy& More table copy$^{\mathrm{a}}$& &  \\
\hline
\multicolumn{4}{l}{$^{\mathrm{a}}$Sample of a Table footnote.}
\end{tabular}
\label{tab1}
\end{center}
\end{table}

\begin{figure}[htbp]
\centerline{\includegraphics{fig1.png}}
\caption{Example of a figure caption.}
\label{fig}
\end{figure}

Figure Labels: Use 8 point Times New Roman for Figure labels. Use words
rather than symbols or abbreviations when writing Figure axis labels to
avoid confusing the reader. As an example, write the quantity
``Magnetization'', or ``Magnetization, M'', not just ``M''. If including
units in the label, present them within parentheses. Do not label axes only
with units. In the example, write ``Magnetization (A/m)'' or ``Magnetization
\{A[m(1)]\}'', not just ``A/m''. Do not label axes with a ratio of
quantities and units. For example, write ``Temperature (K)'', not
``Temperature/K''.

\section*{Acknowledgment}

The preferred spelling of the word ``acknowledgment'' in America is without
an ``e'' after the ``g''. Avoid the stilted expression ``one of us (R. B.
G.) thanks $\ldots$''. Instead, try ``R. B. G. thanks$\ldots$''. Put sponsor
acknowledgments in the unnumbered footnote on the first page.

\section*{References}

Please number citations consecutively within brackets \cite{b1}. The
sentence punctuation follows the bracket \cite{b2}. Refer simply to the reference
number, as in \cite{b3}---do not use ``Ref. \cite{b3}'' or ``reference \cite{b3}'' except at
the beginning of a sentence: ``Reference \cite{b3} was the first $\ldots$''

Number footnotes separately in superscripts. Place the actual footnote at
the bottom of the column in which it was cited. Do not put footnotes in the
abstract or reference list. Use letters for table footnotes.

Unless there are six authors or more give all authors' names; do not use
``et al.''. Papers that have not been published, even if they have been
submitted for publication, should be cited as ``unpublished'' \cite{b4}. Papers
that have been accepted for publication should be cited as ``in press'' \cite{b5}.
Capitalize only the first word in a paper title, except for proper nouns and
element symbols.

For papers published in translation journals, please give the English
citation first, followed by the original foreign-language citation \cite{b6}.


\section{Introduction}\label{intro}
As we know, the linear least squares problem is a classical problem in numerical linear algebra and has numerous important applications in many fields.
There are many direct methods for solving this problem such as
the method by
the QR factorization with pivoting and the method by
the SVD \cite{Osborne1996,Higham2002}. However, these direct methods usually require high storage and are expensive when the involved matrix is large-scale. Hence, some iterative methods are provided for solving large linear least squares problem.
The Gauss-Seidel method \cite{Saad2003} is a very famous one, which,  at each iteration,
first selects a coordinate $j_k\in \{1, \ldots, n\}$ and then minimizes the objective $L(x)=\frac{1}{2} \|{b}-{A} {x}\|^2_{2}$ with respect to this coordinate. This leads to the following iterative process:
\begin{equation*}
\label{gs}
x_{k+1}=x_{k}+\frac{ A^{T}_{(j_k)}(b-Ax_{k}) }{ \| A_{\left(j_{k}\right)} \|_{2}^{2}}e_{j_k},
\end{equation*}
where $A_{\left(j_{k}\right)}$ denotes the $j_k$th column of $A\in R^{m\times n}$ and $e_{j_k}$ is the $j_k$th coordinate basis column vector, i.e., the $j_k$th column of the identity matrix $I$

In 2010,  by using a probability distribution to select the column of $A$ randomly in each iteration,
Leventhal and Lewis \cite{Leventhal2010} considered the randomized Gauss-Seidel (RGS) method, which is also called the randomized coordinate descent method, and showed that it converges linearly in expectation to the solution. Subsequently, many works on RGS method were reported
due to its better performance; see, for example \cite{Ma2015,Hefny2017,Edalatpour2017,Tu2017,Chen2017,Tian2017,Xu1,Razaviyayn2019,Dukui2019} and references therein.

In 2018, Wu \cite{WEIMAGGIE2018CONVERGENCEOT} developed a randomized block Gauss-Seidel (RBGS) method, i.e., a block version of the RGS method,
which
minimizes the objective through multiple coordinates at the same time in each iteration. More specifically, for a subset $\tau \subset \{1, \ldots, n\}$ with $|\tau|=T$, i.e., the number of the elements of the set $\tau$ is $T$, denote the submatrix of $A$ whose columns are indexed by $\tau$ by $A_{\tau}$ and generate a submatrix $I_{\tau}$ of the identity matrix $I$ similarly.
The iterative process can be written as
\begin{equation}
\label{rbgs}
x_{k+1}=x_{k}+I_{\tau}A^{\dagger}_{\tau}(b-Ax_{k}),
\end{equation}
where $A^{\dagger}_{\tau}$ denotes the Moore-Penrose pseudoinverse of $A_{\tau}$. The discussions in \cite{WEIMAGGIE2018CONVERGENCEOT} show that the RBGS method has better convergence compared with the RGS method. Later, Du and Sun \cite{du2019doubly} generalized the RGS method to a doubly stochastic block Gauss-Seidel method.

Recently, by introducing a more efficient probability criterion for selecting the working column from the matrix $A$, Bai and Wu \cite{Bai2019} proposed a greedy randomized coordinate descent (GRCD) method,
which is faster than the RGS method in terms of the number of iterations and computing time.
The idea of greed applied in \cite{Bai2019}  has wide applications and has been used in many works, see for example \cite{Griebel2012,Nguyen2017,Bai2018,Bai2018r,Nutini2018,Zhang2019,Du2019,Liu2019,Haddock2019,Rebrova2019,Niu2020} and references therein.

Inspired by the ideas of the RBGS and GRCD methods,
we develop a
greedy block Gauss-Seidel (GBGS) method in the present paper, and show that
the convergence  factor of the GBGS method can be much smaller than that of the GRCD method.
Experimental results also demonstrate that the GBGS method can significantly accelerate the GRCD method. A drawback of the GBGS method is that, in the block update rule like \cref{rbgs}, we have to
compute the Moore-Penrose pseudoinverses in each iteration, which is quite expensive and also makes the method
be not adequate for distributed implementations. To tackle this problem, inspired by the idea in \cite{Necoara2019,Du20202},
we present a pseudoinverse-free greedy block Gauss-Seidel (PGBGS) method, which
improves the GBGS method in computation cost.

The rest of this paper is organized as follows. In \cref{sec2}, additional notation and preliminaries are provided. We present our GBGS and PGBGS methods and their corresponding convergence properties in \cref{sec3} and \cref{sec4}, respectively. Experimental results are given in \cref{sec5}.
\section{Notation and preliminaries}\label{sec2}

For a vector $z\in R^{n}$, $z^{(j)}$ represents its $j$th entry. For a matrix $G=(g_{ij})\in R^{m\times n}$,
$\sigma_{\min}(G)$, $\sigma_{\max}(G)$, $\|G\|_2$, $\|G\|_F$ and $\mathrm{range}(G)$ denote its
minimum positive singular value, maximum singular value, spectral norm, Frobenius norm, and column space, respectively.
Moreover, if the matrix $G\in R^{n\times n}$ is positive definite, then we define the energy norm of any vector $x\in R^{n}$ as $\| x\|_G:=\sqrt{x^TGx}$.

In the following,
we use $x_{\star}=A^{\dag}b$
to denote the unique least squares solution to the linear least squares problem:
\begin{equation}
\label{ls}
\min \limits _{\mathbf{x} \in \mathbb{R}^{n}}\|{b}-{A} {x}\|^2_{2},
\end{equation}
where $A\in R^{m\times n}$ is of full column rank and $b\in R^{m}$.
It is known that the solution $x_\star$
is the solution to the following normal equation \cite{Osborne1961} for \cref{ls}:
\begin{equation}
\label{1}
A^TAx=A^Tb.
\end{equation}
On the basis of  \cref{1}, Bai and Wu \cite{Bai2019} proposed the GRCD method listed in \cref{alg1}, where $r_k = b-Ax_k$ denotes the residual vector.
\begin{algorithm}
\caption{The GRCD method}
\label{alg1}
\begin{algorithmic}[!htb]
\STATE{\mbox{Input:} ~$A\in R^{m\times n}$, $b\in R^{m}$, $\ell$ , initial estimate $x_0$}
\STATE{\mbox{Output:} ~$x_\ell$}
\FOR{$k=0, 1, 2, \ldots, \ell-1$}
\STATE{Compute
\begin{align}
  \delta_{k}=\frac{1}{2}\left(\frac{1}{\left\|A^Tr_k\right\|_{2}^{2}} \max \limits _{1 \leq j  \leq n}\left\{\frac{\left|A^T_{(j)}r_k\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}}\right\}+\frac{1 }{\|A\|_{F}^{2}}\right).\notag
\end{align}}
\STATE{Determine the index set of positive integers
\begin{align}
 \mathcal{V}_{k}=\left\{j \Bigg|  \left|A^T_{(j)}r_k\right|^{2} \geq \delta_{k} \left\|A^Tr_k\right\|_{2}^{2} \left\|A_{\left(j\right)}\right\|_{2}^{2}\right\}.\notag
\end{align}}
\STATE{Let $s_k=A^Tr_k$ and define $\tilde{s}_k$ as follows \[
\tilde{s}_{k}^{(j)}=\left\{\begin{array}{ll}{s_{k}^{(j)},} & {\text { if } j \in \mathcal{V}_{k}}, \\ {0,} & {\text { otherwise. }}\end{array}\right.
\]}
\STATE{Select $j_k \in \mathcal{V}_{k}$ with probability $Pr(\mathrm{column} = j_k)=\frac{|\tilde{s}_{k}^{(j_k)}|^2}{\|\tilde{s}_{k}\|^2_2}$.}
\STATE{Set
\[x_{k+1}=x_{k}+\frac{ A^{T}_{(j_k)}r_{k} }{ \| A_{\left(j_{k}\right)} \|_{2}^{2}}e_{j_k}.\]}
\ENDFOR
\end{algorithmic}
\end{algorithm}

In addition, to analyze the convergence of our new methods, we need  the following simple result, which can be found in \cite{Bai2018}.
\begin{lemma}
\label{lem1}
For any vector $x \in  \mathrm{range} (A) $, it holds that
\begin{align*}
\|A^{T}x\|^2_2\geq\sigma^{2}_{\min}\left( A\right)\|x\|^2_2.  \label{1210}
\end{align*}
\end{lemma}

\section{Greedy block Gauss-Seidel method}\label{sec3}

There are two main steps in the GBGS method. The first step is to devise a greedy rule to decide the index set $\mathcal{J}_{k}$ whose specific definition is given in \cref{alg2}, and the second step is to update $x_k$ using a update formula like \cref{rbgs}, with which we can
minimize the objective function through all coordinates from the index set $\mathcal{J}_{k}$ at the same time. Note that the GRCD method only updates one coordinate in each iteration.

Based on the above introduction, we propose
\cref{alg2}.
\begin{algorithm}[!htb]
\caption{The GBGS method}
\label{alg2}
\begin{algorithmic}
\STATE{\mbox{Input:} ~$A\in R^{m\times n}$, $b\in R^{m}$, $\theta \in[0,1]$, $\ell$,
 initial estimate $x_0\in R^{n}$}
\STATE{\mbox{Output:} ~$x_\ell$}
\FOR{$k=0, 1, 2, \ldots, \ell-1$}
\STATE{Compute
\begin{align}
\varepsilon_{k}=\frac{\theta}{\left\|A^Tr_k\right\|_{2}^{2}} \max \limits _{1 \leq j  \leq n}\left\{\frac{\left|A^T_{(j)}r_k\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}}\right\}+\frac{1-\theta }{\|A\|_{F}^{2}} .\notag
\end{align}}
\STATE{Determine the index set of positive integers
\begin{align}
 \mathcal{J}_{k}=\left\{j_k \Bigg|  \left|A^T_{(j_k)}r_k\right|^{2} \geq \varepsilon_{k}\left\|A^Tr_k\right\|_{2}^{2}  \left\|A_{\left(j_k\right)}\right\|_{2}^{2}\right\}.\notag
\end{align}}
\STATE{Set
\begin{align}
\label{up2}
x_{k+1}=x_{k}+I_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}}r_{k}.
\end{align}}
\ENDFOR
\end{algorithmic}
\end{algorithm}

\begin{remark}
\label{rmk1}
As done in \cite{Bai2018r}, if we replace $ \delta_{k}$ in \cref{alg1} by $\varepsilon_{k}$ in \cref{alg2}, we can get a relaxed version of the GRCD method.
In addition, it is easy to see that if $\theta=\frac{1}{2}$, then $\mathcal{J}_{k}=\mathcal{V}_{k}$, i.e., the index sets of the GBGS and GRCD methods are the same.
\end{remark}

\begin{remark}
\label{rmk2}
Note that if
\[\frac{\left|A^T_{(j_k)}r_k\right|^{2}}{\left\|A_{\left(j_k\right)}\right\|_{2}^{2}}=  \max \limits _{1 \leq j  \leq n}\left\{\frac{\left|A^T_{(j)}r_k\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}}\right\},\]
then $j_k\in\mathcal{J}_{k}.$ So the index set $\mathcal{J}_{k}$ in \cref{alg2} is nonempty for all iteration index $k$.
\end{remark}

In the following, we give the convergence theorem of the GBGS method.

\begin{theorem}
\label{theorem1}
The iteration sequence $\{x_k\}_{k=0}^\infty$ generated by \cref{alg2}, starting from an initial guess $x_0\in R^{n}$, converges linearly to the unique least squares solution $x_{\star}=A^{\dag}b$ and
\begin{equation}
\label{4}
  \| x_1-x_\star\|^2_{A^TA} \leq\left(1- \frac{\|A_{\mathcal{J}_{0}}\|^2_F}{\sigma^2_{\max}(A_{\mathcal{J}_{0}})}  \frac{\sigma^{2}_{\min}(A)}{\left\|A\right\|_{F}^{2}} \right) \| x_0-x_\star\|^2_{A^{T}A},
\end{equation}
and
 \begin{align}
\label{5}
\| x_{k+1}-x_\star\|^2_{A^TA} \leq\left(1- \frac{\|A_{\mathcal{J}_{k}}\|^2_F}{\sigma^2_{\max}(A_{\mathcal{J}_{k}})} \left(\theta\frac{\|A\|^2_F}{\|A\|^2_F-\|A_{\mathcal{J}_{k-1}}\|^2_F}+(1-\theta)\right) \frac{\sigma^{2}_{\min}(A)}{\|A\|^2_F}\right)\| x_k-x_\star\|^2_{A^TA},\\
k=1, 2, \ldots.\notag
\end{align}
Moreover, let $\alpha=\max \{\sigma^2_{\max}(A_{\mathcal{J}_{k}})\}$, $\beta=\min\{\|A_{\mathcal{J}_{k}}\|^2_F\}$, and $\gamma=\max\{\|A\|^2_F-\|A_{\mathcal{J}_{k-1}}\|^2_F\}$.
Then
\begin{align}
\label{6}
\| x_{k}-x_\star\|^2_{A^TA} &\leq\left(1- \frac{\beta}{\alpha}(\theta\frac{\|A\|^2_F}{\gamma}+(1-\theta)) \frac{\sigma^{2}_{\min}(A)}{\|A\|^2_F} \right)^{k-1}\\
&\times\left(1- \frac{\|A_{\mathcal{J}_{0}}\|^2_F}{\sigma^2_{\max}(A_{\mathcal{J}_{0}})}\frac{\sigma^{2}_{\min}(A)}{\left\|A\right\|_{F}^{2}} \right)\| x_0-x_\star\|^2_{A^TA},
\quad k=1, 2, \ldots.\notag
\end{align}
\end{theorem}
\begin{proof}
 From the update formula \cref{up2} in \cref{alg2}, we have
 \begin{align}
A(x_{k+1}-x_{\star})&=~A(x_{k}-x_{\star})+A_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}}r_{k},  \notag
  \end{align}
  which together with the fact $A_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}}=(A_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}})^{T}$ and $A^TAx_{\star}=A^Tb$ gives
\begin{align}
A(x_{k+1}-x_{\star})&= ~ A(x_{k}-x_{\star})+(A^{\dag}_{\mathcal{J}_{k}})^{T}A^{T}_{\mathcal{J}_{k}}(b-Ax_{k})  \notag\\
&= ~\left(I-A_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}} \right)A\left(x_{k}-x_\star\right). \notag
\end{align}
Since $A_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}}$ is an orthogonal projector, taking the square of the Euclidean norm on both sides of the above equation and applying Pythagorean theorem, we get
\begin{align}
\|A(x_{k+1}-x_{\star})\|^2_2&= ~\|\left(I-A_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}} \right)A\left(x_{k}-x_\star\right)\|^2_2\notag\\
&= ~\|A\left(x_{k}-x_\star\right)\|^2_2-\|A_{\mathcal{J}_{k}}A^{\dag}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2\notag\\
&= ~\|A\left(x_{k}-x_\star\right)\|^2_2-\|(A^{\dag}_{\mathcal{J}_{k}})^{T}A^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2,\notag
\end{align}
or equivalently,
\begin{align}
\left\|x_{k+1}-x_{\star}\right\|^2_{A^{T}A}=\left\|x_{k}-x_{\star}\right\|^2_{A^{T}A}-\|(A^{\dag}_{\mathcal{J}_{k}})^{T}A^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2,\notag
\end{align}
which together with \cref{lem1} and the fact $\sigma^2_{\min}(A^{\dag}_{\mathcal{J}_{k}})=\sigma^{-2}_{\max}(A_{\mathcal{J}_{k}})$ yields
\begin{align}
\left\|x_{k+1}-x_{\star}\right\|^2_{A^{T}A}&\leq~\left\|x_{k}-x_{\star}\right\|^2_{A^{T}A}-\sigma^2_{\min}(A^{\dag}_{\mathcal{J}_{k}})\|A^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2\notag\\
&=~\left\|x_{k}-x_{\star}\right\|^2_{A^{T}A}-\frac{1}{\sigma^2_{\max}(A_{\mathcal{J}_{k}})}\|A^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2.\notag
\end{align}
On the other hand, from \cref{alg2}, we have
\begin{align}
\|A^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2= \sum\limits_{j_k\in  \mathcal{J}_{k}}  \left|A^T_{(j_k)}r_k\right|^2
\geq\sum\limits_{j_k\in  \mathcal{J}_{k}} \varepsilon_{k} \left\|A^Tr_k\right\|_{2}^{2} \left\|A_{\left(j_k\right)}\right\|_{2}^{2}
= \varepsilon_{k}\left\|A^Tr_k\right\|_{2}^{2}  \|A_{\mathcal{J}_{k}}\|^2_F. \notag
\end{align}
Then,
\begin{align}
\left\|x_{k+1}-x_{\star}\right\|^2_{A^{T}A}\leq\left\|x_{k}-x_{\star}\right\|^2_{A^{T}A}-\frac{\|A_{\mathcal{J}_{k}}\|^2_F}{\sigma^2_{\max}(A_{\mathcal{J}_{k}})}\left\|A^Tr_k\right\|_{2}^{2} \varepsilon_{k}.\label{1217}
\end{align}

For $k=0$, we have
\begin{align}
\left\|A^Tr_0\right\|_{2}^{2} \varepsilon_{0}&=~\theta \max \limits _{1 \leq j  \leq n}\left\{\frac{\left|A^T_{(j)}r_0\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}}\right\}+(1-\theta )\frac{ \left\|A^Tr_0\right\|_{2}^{2}}{\|A\|_{F}^{2}}\notag
\\
&\geq~\theta \sum\limits_{j=1}^{n} \frac{\left\|A_{\left(j\right)}\right\|_{2}^{2}}{\|A\|_{F}^{2}} \frac{\left|A^T_{(j)}r_0\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}}+(1-\theta )\frac{ \left\|A^Tr_0\right\|_{2}^{2}}{\|A\|_{F}^{2}}\notag \\
&=~\theta\frac{\|A^Tr_0\|^{2}_2}{\left\|A\right\|_{F}^{2}}+(1-\theta )\frac{ \left\|A^Tr_0\right\|_{2}^{2}}{\|A\|_{F}^{2}}=\frac{ \left\|A^Tr_0\right\|_{2}^{2}}{\|A\|_{F}^{2}},\notag
\end{align}
which together with \cref{lem1} leads to
\begin{align}
\left\|A^Tr_0\right\|_{2}^{2} \varepsilon_{0}\geq  \frac{1}{\left\|A\right\|_{F}^{2}} \sigma^{2}_{\min}(A)\|A(x_0-x_\star)\|^{2}_2=  \frac{\sigma^{2}_{\min}(A)}{\left\|A\right\|_{F}^{2}} \|x_0-x_\star\|^{2}_{A^{T}A}.\label{1212111}
\end{align}
Thus, substituting \cref{1212111} into \cref{1217}, we obtain
\begin{align}
\left\|x_{1}-x_{\star}\right\|^2_{A^{T}A}&\leq~\left\|x_{0}-x_{\star}\right\|^2_{A^{T}A}-\frac{\|A_{\mathcal{J}_{0}}\|^2_F}{\sigma^2_{\max}(A_{\mathcal{J}_{0}})}  \frac{\sigma^{2}_{\min}(A)}{\left\|A\right\|_{F}^{2}} \|x_0-x_\star\|^{2}_{A^{T}A} \notag\\
&=~\left(1- \frac{\|A_{\mathcal{J}_{0}}\|^2_F}{\sigma^2_{\max}(A_{\mathcal{J}_{0}})}  \frac{\sigma^{2}_{\min}(A)}{\left\|A\right\|_{F}^{2}} \right) \| x_0-x_\star\|^2_{A^{T}A}, \notag
\end{align}
which is just the estimate \cref{4}.

For $k\geq1$, to find the lower bound of $\left\|A^Tr_k\right\|_{2}^{2} \varepsilon_{k}$, first note that
\begin{align}
 A^T_{ \mathcal{J}_{k-1} }r_k&=  ~A^T_{\mathcal{J}_{k-1}}\left(b-A(x_{k-1}+I_{\mathcal{J}_{k-1}}A^{\dag}_{\mathcal{J}_{k-1}}r_{k-1})  \right)   \notag\\
&= ~A^T_{ \mathcal{J}_{k-1}}r_{k-1}-A^T_{ \mathcal{J}_{k-1}} A_{\mathcal{J}_{k-1}}A^{\dag}_{\mathcal{J}_{k-1}}r_{k-1}  \notag\\
&= ~0, \notag
\end{align}
where we have used the update formula \cref{up2} and the property of the Moore-Penrose pseudoinverse.
Then
\begin{align}
\|A^{T}r_k\|^2_2&=~\sum\limits_{j=1}^{n}|A^{T}_{(j)}r_k|^2=\sum\limits_{  j=1 \atop  j\notin  \mathcal{J}_{k-1}}^n|A^{T}_{(j)}r_k|^2 =~\sum\limits_{  j=1 \atop  j\notin  \mathcal{J}_{k-1}}^n\frac{|A^{T}_{(j)}r_k|^2}{\|A_{j}\|^2_2} \|A_{j}\|^2_2\notag\\
&\leq~\max\limits_{1\leq j\leq n  }\frac{|A^{T}_{(j)}r_k|^2}{\|A_{j}\|^2_2}  \sum\limits_{  j=1 \atop  j\notin  \mathcal{J}_{k-1}}^n\|A_{j}\|^2_2=~\max\limits_{1\leq j\leq n  }\frac{|A^{T}_{(j)}r_k|^2}{\|A_{j}\|^2_2} (\|A\|^2_F-\|A_{\mathcal{J}_{k-1}}\|^2_F),\notag
\end{align}
which can first imply a lower bound of $\max\limits_{1\leq j\leq n  }\frac{|A^{T}_{(j)}r_k|^2}{\|A_{j}\|^2_2}$ and then of $\left\|A^Tr_k\right\|_{2}^{2} \varepsilon_{k}$ as follows
\begin{align}
\left\|A^Tr_k\right\|_{2}^{2} \varepsilon_{k
\geq~ \theta\frac{\|A^{T}r_k\|^2_2}{\|A\|^2_F-\|A_{\mathcal{J}_{k-1}}\|^2_F}+(1-\theta)\frac{ \left\|A^Tr_k\right\|_{2}^{2}}{\|A\|_{F}^{2}}.\notag
\end{align}
Further, considering
\cref{lem1}, we have
\begin{align}
\left\|A^Tr_k\right\|_{2}^{2}\varepsilon_{k}&\geq~  \left(\theta\frac{1}{\|A\|^2_F-\|A_{\mathcal{J}_{k-1}}\|^2_F}+(1-\theta)\frac{ 1}{\|A\|_{F}^{2}}\right) \sigma^{2}_{\min}(A) \|A(x_k-x_\star)\|^2_2\notag\\
&=~ \left(\theta\frac{\|A\|^2_F}{\|A\|^2_F-\|A_{\mathcal{J}_{k-1}}\|^2_F}+(1-\theta)\right) \frac{\sigma^{2}_{\min}(A)}{\|A\|^2_F} \|x_k-x_\star\|^{2}_{A^{T}A}.\label{12121111}
\end{align}
Thus, substituting \cref{12121111} into \cref{1217} gives
the estimate \cref{5}. By induction on the iteration index $k$, we can obtain the estimate \cref{6}.
\end{proof}

\begin{remark}
\label{rmk3111}
If $\theta=\frac{1}{2}$, i.e., the index sets of the GBGS and GRCD methods are the same, the first term in the right side of \cref{5} reduces to
$$
\eta=1- \frac{\|A_{\mathcal{J}_{k}}\|^2_F}{\sigma^2_{\max}(A_{\mathcal{J}_{k}})} \frac{1}{2}\left( \frac{\|A\|^2_F}{\|A\|^2_F-\|A_{\mathcal{J}_{k-1}}\|^2_F}+1\right) \frac{\sigma^{2}_{\min}(A)}{\|A\|^2_F}.
$$
Since
\[ \|A\|^2_{F}-\|A_{\mathcal{J}_{k-1}}\|^2_{F}\leq\|A\|^2_F-\min\limits_{1\leq j\leq n}\|A_{(j)}\|^{2}_2\ \textrm{ and }\ \frac{\|A_{\mathcal{J}_{k}}\|^2_{F}}{\sigma^2_{\max}(A_{\mathcal{J}_{k}})}\geq1,\]
we have
\begin{eqnarray*}
\eta \leqslant 1- \frac{1}{2} \left(\frac{\|A\|_{F}^{2}}{\|A\|_{F}^{2}-\min \limits_{1 \leq j \leq n}\left\|A_{(j)}\right\|_{2}^{2}}+1 \right)\frac{\sigma^2_{\min}(A)}{\|A\|_{F}^{2}}.
\end{eqnarray*}
Note that the error estimate in expectation of the GRCD method given in \cite{Bai2019} is
\[
\mathbb{E}_{k}\left\|x_{k+1}-x_{\star}\right\|_{A^{T} A}^{2} \leq\left (1- \frac{1}{2} \left(\frac{\|A\|_{F}^{2}}{\|A\|_{F}^{2}-\min \limits_{1 \leq j \leq n}\left\|A_{(j)}\right\|_{2}^{2}}+1 \right)\frac{\sigma^2_{\min}(A)}{\|A\|_{F}^{2}}\right  )\left\|x_{k}-x_{\star}\right\|_{A^{T} A}^{2} , \]
where $ k=1,2, \ldots.$ So the convergence  factor of the GBGS method is
smaller than that of the GRCD method in this case,
and the former can be much smaller than the latter because
$\|A\|^2_{F}-\|A_{\mathcal{J}_{k-1}}\|^2_{F}$ can be much smaller than $\|A\|^2_F-\min\limits_{1\leq j\leq n}\|A_{(j)}\|^{2}_2$ and $\|A_{\mathcal{J}_{k}}\|^2_{F}$ can be much larger than $\sigma^2_{\max}(A_{\mathcal{J}_{k}})$.

\end{remark}

\begin{remark}
\label{rmk4}
Very recently, Niu and Zheng \cite{Niu2020} provided a greedy block Kaczmarz algorithm for solving large-scale linear systems, which is a block version of the famous greedy randomized Kaczmarz method given in \cite{Bai2018}. The technique in \cite{Niu2020} can be also applied to develop a block version of the greedy Gauss-Seidel method.

\end{remark}

\section{Pseudoinverse-free greedy block Gauss-Seidel method}\label{sec4}
As mentioned in \cref{intro},  there is a drawback in the update rule \cref{up2} in \cref{alg2}, that is,
we need to compute the Moore-Penrose pseudoinverse of the column submatrix $A_{\mathcal{J}_{k}}$ in each iteration.
Inspired by the pseudoinverse-free block Kaczmarz methods \cite{Necoara2019, Du20202},
we design a PGBGS method.
More specifically, assume
 $|\mathcal{J}_{k}|=t$ and let $\mathcal{J}_{k}=\{j_{k1}, j_{k2}, \ldots, j_{kt}\}$.
 Rather than using the update rule \cref{up2}, we execute all coordinates in $\mathcal{J}_{k}$ simultaneously using the following update rule
\begin{align}
x_{k+1}&=~x_{k}+\omega\frac{ A^{T}_{(j_{k1})}r_{k} }{ \| A_{\left(j_{k1}\right)} \|_{2}^{2}}e_{j_{k1}}+\omega\frac{ A^{T}_{(j_{k2})}r_{k} }{ \| A_{\left(j_{k2}\right)} \|_{2}^{2}}e_{j_{k2}}+\cdots+\omega\frac{ A^{T}_{(j_{kt})}r_{k} }{ \| A_{\left(j_{kt}\right)} \|_{2}^{2}}e_{j_{kt}} \notag\\
&=~x_{k}+\omega\frac{e_{j_{k1}}}{\| A_{\left(j_{k1}\right)} \|_{2}}\frac{ A^{T}_{(j_{k1})}}{ \| A_{\left(j_{k1}\right)} \|_{2}}r_{k}+\omega\frac{e_{j_{k2}}}{\| A_{\left(j_{k2}\right)} \|_{2}}\frac{ A^{T}_{(j_{k2})}}{ \| A_{\left(j_{k2}\right)} \|_{2}}r_{k}+\cdots+\omega\frac{e_{j_{kt}}}{\| A_{\left(j_{kt}\right)} \|_{2}}\frac{ A^{T}_{(j_{kt})}}{ \| A_{\left(j_{kt}\right)} \|_{2}}r_{k}\notag\\
&=~x_{k}+\omega\widetilde{I}_{\mathcal{J}_{k}}\widetilde{A}^{T}_{\mathcal{J}_{k}}r_{k},\label{derive}
\end{align}
where
\begin{eqnarray*}
\widetilde{I}_{\mathcal{J}_{k}}&=&\left[\frac{e_{j_{k1}}}{\| A_{\left(j_{k1}\right)} \|_{2}}, \frac{e_{j_{k2}}}{\| A_{\left(j_{k2}\right)} \|_{2}}, \ldots, \frac{e_{j_{kt}}}{\| A_{\left(j_{kt}\right)} \|_{2}}\right]\in R^{n\times t},\\
\widetilde{A}^{T}_{\mathcal{J}_{k}}&=&\left[\frac{ A_{(j_{k1})}}{ \| A_{\left(j_{k1}\right)} \|_{2}}, \frac{ A_{(j_{k2})}}{ \| A_{\left(j_{k2}\right)} \|_{2}}, \ldots, \frac{ A_{(j_{kt})}}{ \| A_{\left(j_{kt}\right)} \|_{2}}\right]^{T}\in R^{t\times m},
\end{eqnarray*}
and $\omega$  is used to adjust the convergence of the method and satisfies
\[
\left\{\begin{array}{ll}{0<\omega<\frac{2}{\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}})},} & {\text { if } \Delta<0}, \\ {0<\omega<\frac{1-\sqrt{\Delta}}{\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}})}~ \text {or} ~ \frac{1+\sqrt{\Delta}}{\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}})}<\omega<\frac{2}{\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}})},} & {\text { otherwise, }}\end{array}\right.
\]
where $\Delta=1-\frac{\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}})\|A\|^2_F}{|\mathcal{J}_{k}|\sigma^{2}_{\min}(A)}$. It should be pointed out that the above conditions are only sufficient  but not necessary ones.

In summary,  we have \cref{alg3}.
\begin{algorithm}[!htb]
\caption{The PGBGS method}
\label{alg3}
\begin{algorithmic}
\STATE{\mbox{Input:} ~$A\in R^{m\times n}$, $b\in R^{m}$, $\theta \in[0,1]$, $\omega$, $\ell$,
initial estimate $x_0\in R^{n}$}
\STATE{\mbox{Output:} ~$x_\ell$}
\FOR{$k=0, 1, 2, \ldots, \ell-1$}
\STATE{Compute
\begin{align}
\varepsilon_{k}=\frac{\theta}{\left\|A^Tr_k\right\|_{2}^{2}} \max \limits _{1 \leq j  \leq n}\left\{\frac{\left|A^T_{(j)}r_k\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}}\right\}+\frac{1-\theta }{\|A\|_{F}^{2}} .\notag
\end{align}}
\STATE{Determine the index set of positive integers
\begin{align}
 \mathcal{J}_{k}=\left\{j_k \Bigg|  \left|A^T_{(j_k)}r_k\right|^{2} \geq \varepsilon_{k}\left\|A^Tr_k\right\|_{2}^{2}  \left\|A_{\left(j_k\right)}\right\|_{2}^{2}\right\}.\notag
\end{align}}
\STATE{Set
\begin{align}
\label{up3}
x_{k+1}=x_{k}+\omega\widetilde{I}_{\mathcal{J}_{k}}\widetilde{A}^{T}_{\mathcal{J}_{k}}r_{k}.
\end{align}}
\ENDFOR
\end{algorithmic}
\end{algorithm}

\begin{remark}
\label{rmk5}
To compare the computation cost of the GRCD, GBGS and PGBGS methods, we present the flops of the update rules of the three methods in \cref{TAB1}, where we have used the fact $A^{\dag}_{\mathcal{J}_{k}}=(A^T_{\mathcal{J}_{k}}A_{\mathcal{J}_{k}})^{-1}A^T_{\mathcal{J}_{k}}$.

\begin{table}[!htbp]\centerin
\begin{small}\scriptsize
\caption{Flops for the update rules of the \texttt{GRCD}, \texttt{GBGS} and \texttt{PGBGS} methods.}\centering \label{TAB1}
 \begin{tabular}{|c|c|c|c|}
 \hline
\textbf{Method}  & \textbf{GRCD}     & \textbf{GBGS} &\textbf{PGBGS} \\
\hline
\textbf{Flops}  & $2(m+n)$     &$2|\mathcal{J}_{k}|^3+(4m-3)|\mathcal{J}_{k}|^2+(m+2n)|\mathcal{J}_{k}|$   &$(2m+2n+1)|\mathcal{J}_{k}|$ \\
\hline
\end{tabular}
\end{small}
\end{table}

It is easy to find that $2|\mathcal{J}_{k}|^3+(4m-3)|\mathcal{J}_{k}|^2+(m+2n)|\mathcal{J}_{k}|>(2m+2n+1)|\mathcal{J}_{k}|$.
So the GBGS method needs more flops compared with the PGBGS method.
Thus, for the same accuracy, the latter needs less computing time. More importantly, from the derivation of \cref{derive}, we can find that the PGBGS  method can be used for distributed implementations, which can yield more significant improvements in computation cost.
Although the GRCD method requires fewer flops in each iteration compared with our two block methods, the GRCD method needs much more iterations because it only executes one column in one iteration. So our methods behave better in the total computing time. Moreover, if $|\mathcal{J}_{k}|$ is small, the difference between the flops needed in the GRCD and PGBGS may be negligible.  These results are confirmed by experimental results given in \cref{sec5}.
\end{remark}

In the following, we give the convergence theorem of the PGBGS method.
\begin{theorem}
\label{theorem2}
The iteration sequence $\{x_k\}_{k=0}^\infty$ generated by \cref{alg3}, starting from an initial guess $x_0\in R^{n}$, converges linearly to the unique least squares solution $x_{\star}=A^{\dag}b$ and
\begin{align}
\left\|x_{k+1}-x_{\star}\right\|^2_{A^{T}A}&\leq~\left(1-(2\omega-\omega^2\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}})) |\mathcal{J}_{k}| \frac{\sigma^{2}_{\min}(A)}{\left\|A\right\|_{F}^{2}}\right)\|x_k-x_\star\|^{2}_{A^{T}A}.\label{4sd}
\end{align}
\end{theorem}

\begin{proof}
From the update rule \cref{up3} in \cref{alg3}, we have
\begin{align}
A(x_{k+1}-x_{\star})&=~A(x_{k}-x_{\star})+\omega\widetilde{A}_{\mathcal{J}_{k}}\widetilde{A}^{T}_{\mathcal{J}_{k}}r_{k},  \notag
\end{align}
which together with the fact $A^TAx_{\star}=A^Tb$ gives
\begin{align}
A(x_{k+1}-x_{\star})&=~\left(I-\omega\widetilde{A}_{\mathcal{J}_{k}}\widetilde{A}^{T}_{\mathcal{J}_{k}} \right)A\left(x_{k}-x_\star\right). \notag
\end{align}
Taking the square of the Euclidean norm on both sides, we get
\begin{align}
\|A(x_{k+1}-x_{\star})\|^2_2&= ~\|\left(I-\omega\widetilde{A}_{\mathcal{J}_{k}}\widetilde{A}^{T}_{\mathcal{J}_{k}} \right)A\left(x_{k}-x_\star\right)\|^2_2\notag\\
&= ~\|A\left(x_{k}-x_\star\right)\|^2_2-2\omega\|\widetilde{A}^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2+\omega^2\|\widetilde{A}_{\mathcal{J}_{k}}\widetilde{A}^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2\notag\\
&\leq~\|A\left(x_{k}-x_\star\right)\|^2_2-(2\omega-\omega^2\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}}))\|\widetilde{A}^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2,\notag
\end{align}
or equivalently,
\begin{equation}
\label{754r}
\left\|x_{k+1}-x_{\star}\right\|^2_{A^{T}A}\leq\left\|x_{k}-x_{\star}\right\|^2_{A^{T}A}-(2\omega-\omega^2\sigma^2_{\max}(\widetilde{A}_{\mathcal{J}_{k}}))\|\widetilde{A}^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2.
\end{equation}
On the other hand, from \cref{alg3} and \cref{derive}, we have
\begin{align}
\|\widetilde{A}^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2
&=~\sum\limits_{j_k\in  \mathcal{J}_{k}} \left|\widetilde{A}^T_{(j_k)}r_k\right|^2=\sum\limits_{j_k\in  \mathcal{J}_{k}} \left|\frac{A^T_{(j_k)}}{\|A_{(j_k)\|_2}}r_k\right|^2\notag\\
&\geq~\sum\limits_{j_k\in  \mathcal{J}_{k}}\frac{1}{\left\|A_{\left(j_k\right)}\right\|_{2}^{2}} \varepsilon_{k} \left\|A^Tr_k\right\|_{2}^{2} \left\|A_{\left(j_k\right)}\right\|_{2}^{2}\notag\\
&=~\sum\limits_{j_k\in  \mathcal{J}_{k}} \left(  \theta  \max \limits _{1 \leq j  \leq n}\left\{\frac{\left|A^T_{(j)}r_k\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}}\right\}+(1-\theta)\frac{ \left\|A^Tr_k\right\|_{2}^{2}}{\|A\|_{F}^{2}} \right) \notag\\
&\geq~|\mathcal{J}_{k}|\left(\theta\sum\limits_{j=1}^n \frac{\left|A^T_{(j)}r_k\right|^{2}}{\left\|A_{\left(j\right)}\right\|_{2}^{2}} \frac{\left\|A_{\left(j\right)}\right\|_{2}^{2}}{\|A\|^2_F}+(1-\theta)\frac{ \left\|A^Tr_k\right\|_{2}^{2}}{\|A\|_{F}^{2}}\right)\notag\\
&=~|\mathcal{J}_{k}|\frac{\|A^Tr_k\|^{2}_{2}}{\|A\|^2_F},\notag
\end{align}
which together with \cref{lem1} gives
\begin{align}
\|\widetilde{A}^{T}_{\mathcal{J}_{k}}A\left(x_{k}-x_\star\right)\|^2_2&\geq~|\mathcal{J}_{k}|   \frac{\sigma^{2}_{\min}(A)}{\left\|A\right\|_{F}^{2}} \|x_k-x_\star\|^{2}_{A^{T}A}.\label{12ww111}
\end{align}
Thus, substituting \cref{12ww111} into \cref{754r}, we obtain
the estimate \cref{4sd}.
\end{proof}

\section{Experimental results}\label{sec5}
In this section, we compare the GRCD, GBGS, and PGBGS methods
using the matrix $A\in R^{m\times n}$ from two sets. One is generated randomly by using the MATLAB function \texttt{randn}, and the other one contains the matrices \texttt{abtaha1} and \texttt{ash958} from the University of Florida sparse matrix collection \cite{Davis2011}. To compare these methods fairly,
we set $\theta=\frac{1}{2}$ in $\varepsilon_{k}$, i.e., let the index sets of the three methods be the same. In addition, we set $\omega=1$ in the PGBGS method.

We compare the three methods mainly in terms of the flops
(denoted as ``Flops''), the iteration numbers (denoted as ``Iteration'') and the computing time in seconds (denoted as ``CPU time(s)''). In our specific experiments, the solution vector $x_\star$ is generated randomly by the MATLAB function \texttt{randn}. For the consistent problem, we set the right-hand side $b=Ax_{\star}$. For the inconsistent problem, we set the right-hand side $b=Ax_{\star}+r_{0}$, where $r_0$ is a nonzero vector belonging to the null space of $A^{T}$ generated by the MATLAB function \texttt{null}. All  the test problems are started from an initial zero vector $x_{0}=0$ and terminated once the \emph{relative solution error} (\texttt{RES}), defined by \[\texttt{RES}=\frac{\left\|x_{k}-x_{\star}\right\|^{2}_2}{\left\|x_{\star}\right\|^{2}_2},\] satisfies $\texttt{RES}<10^{-6}$ or the number of iteration steps exceeds $ 200,000$.

For the first class of matrices, that is, the matrices generated randomly, the numerical results on RES in base-10 logarithm versus the Flops, Iteration and CPU time(s) of the three methods are presented in \cref{fig1} when the linear system is consistent, and in \cref{fig2} when the linear system is inconsistent.

\begin{figure}[!htb]
 \begin{center}
\includegraphics [height=4cm,width=15cm ]{consistent_5000_1000.eps}
\includegraphics [height=4cm,width=15cm ]{consistent_5000_2000.eps}
 \end{center}
\caption{\rm Results for dense consistent linear systems with matrices generated randomly. (up) $A$ is of order $5000\times 1000$; (down) $A$ is of order $5000\times 2000$.}\label{fig1}
\end{figure}
\begin{figure}[!htb]
 \begin{center}
\includegraphics [height=4cm,width=15cm ]{inconsistent_5000_1000.eps}
\includegraphics [height=4cm,width=15cm ]{inconsistent_5000_2000.eps}
 \end{center}
\caption{\rm Results for dense inconsistent linear systems with matrices generated randomly. (up) $A$ is of order $5000\times 1000$; (down) $A$ is of order $5000\times 2000$.}\label{fig2}
\end{figure}

From \cref{fig1,fig2}, we can find that the PGBGS method requires almost the same number of flops as the GRCD method, and the GBGS method needs the most flops.
However, the GBGS method requires the fewest iterations, and, as desired, the PGBGS method needs the least computing time. Moreover, the differences in iterations and computing time between our methods and the GRCD method are remarkable. These results are consistent with the analysis in \cref{rmk5}.

For the second class of matrices, that is, the sparse full column rank matrices from \cite{Davis2011}, we plot the numerical results on RES in base-10 logarithm versus the Flops, Iteration and CPU time(s) of the three methods in \cref{fig3} when the linear system is consistent, and in \cref{fig4} when the linear system is inconsistent.

\begin{figure}[!htb]
 \begin{center}
\includegraphics [height=4cm,width=15cm ]{consistent_abtaha1.eps}
\includegraphics [height=4cm,width=15cm ]{consistent_ash958.eps}
 \end{center}
\caption{\rm Results for sparse consistent linear systems with matrices from the University of Florida sparse matrix collection. (up) \texttt{abtaha1}; (down) \texttt{ash958} .}\label{fig3}
\end{figure}
\begin{figure}[!htb]
 \begin{center}
\includegraphics [height=4cm,width=15cm ]{inconsistent_abtaha1.eps}
\includegraphics [height=4cm,width=15cm ]{inconsistent_ash958.eps}
 \end{center}
\caption{\rm Results for sparse inconsistent linear systems  with matrices from the University of Florida sparse matrix collection. (up) \texttt{abtaha1}; (down) \texttt{ash958}.}\label{fig4}
\end{figure}

The findings from \cref{fig3,fig4} are similar to the ones from \cref{fig1,fig2}. That is, the GBGS and PGBGS methods can far outperform the GRCD method in terms of the number of iterations and computing time.

Therefore, in all the cases, although the GBGS method needs more flops because of the computation of the Moore-Penrose pseudoinverse, it is still much faster than the GRCD method. This is mainly because the latter only updates one coordinate in each iteration while the former executes multiple coordinates. The PGBGS method has more advantages because it doesn't need to calculate the Moore-Penrose pseudoinverse any more.

 \clearpage
\bibliographystyle{siamplain}