%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%  [doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

%\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
\documentclass[linenumbers, doublespacing]{bmcart}

%%% Load packages
\usepackage{amsthm,amsmath}
%\RequirePackage{natbib}
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
%\RequirePackage{hyperref}
\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\includegraphic{}
\def\includegraphics{}

%%% Put your definitions there:
\startlocaldefs

\usepackage[pdftex]{graphicx}
%\usepackage{subfig}
\usepackage{amssymb}
%\usepackage{breqn}
%\usepackage{latexsym}
\usepackage[linesnumbered, ruled]{algorithm2e}

\usepackage{array, makecell} 

\usepackage{etoolbox}
\makeatletter
\patchcmd\set@numberlines@box{\llap}{\@gobble}{}{}
\patchcmd\set@numberlines@box
{\rlap{\numberlines@size\kern\numberlines@skip\the\@tempcnta}}
{\rlap{\numberlines@size\kern 10pt\the\@tempcnta}}{}{}
\makeatother
%\usepackage[right, running]{lineno}
%\linenumbers

\def\*#1{\mathbf{#1}}
\DeclareTextFontCommand{\textbfit}{%
  \fontseries\bfdefault % change series without selecting the font yet
  \itshape
}
\makeatletter
\renewcommand*\env@matrix[1][c]{\hskip -\arraycolsep
  \let\@ifnextchar\new@ifnextchar
  \array{*\c@MaxMatrixCols #1}}
\makeatother

\newcommand*{\defeq}{\stackrel{\text{def}}{=}}

\endlocaldefs



%%% Begin ...
\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Orthogonal Approach to Independent Component Analysis Using Quaternionic Factorization}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author[
   addressref={aff1},                   % id's of addresses, e.g. {aff1,aff2}
   corref={aff1},                       % id of corresponding address, if any
%   noteref={n1},                        % id's of article notes, if any
   email={a.borowicz@pb.edu.pl}   % email address
]{\inits{AB}\fnm{Adam} \snm{Borowicz}}
%\author[
%   addressref={aff1,aff2},
%   email={john.RS.Smith@cambridge.co.uk}
%]{\inits{JRS}\fnm{John RS} \snm{Smith}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\address[id=aff1]{%                           % unique id
  \orgname{Faculty of Computer Science, Bialystok University of Technology}, % university, etc
  \street{Wiejska 45A},                     %
  \postcode{15-351}                                % post or zip code
  \city{Bialystok},                              % city
  \cny{Poland}                                    % country
}
%\address[id=aff2]{%
%  \orgname{Marine Ecology Department, Institute of Marine Sciences Kiel},
%  \street{D\"{u}sternbrooker Weg 20},
%  \postcode{24105}
%  \city{Kiel},
%  \cny{Germany}
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{artnotes}
%\note{Sample of title note}     % note to the article
%\note[id=n1]{Equal contributor} % note, connected to author
\end{artnotes}

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract} % abstract
Independent component analysis (ICA) is a popular technique for demixing multi-channel data.
The performance of a typical ICA algorithm strongly depends on the presence of additive noise, the actual distribution of source signals, and the estimated number of non-Gaussian components.
Often a linear mixing model is assumed and source signals are extracted by data whitening followed by a sequence of plane (Jacobi) rotations.
In this article, we develop a novel algorithm, based on the quaternionic factorization of rotation matrices and the Newton-Raphson iterative scheme.
Unlike conventional rotational techniques such as the JADE algorithm, our method exploits $4 \times 4$ rotation matrices and uses approximate negentropy as a contrast function.
Consequently, the proposed method can be adjusted to a given data distribution (e.g. super-Gaussians) by selecting a suitable non-linear function that approximates the negentropy.
Compared to the widely-used, the symmetric FastICA algorithm, the proposed method does not require an orthogonalization step and is more accurate in the presence of multiple Gaussian sources.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The keywords begin here                  %%
%%                                          %%
%% Put each keyword in separate \kwd{}.     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{keyword}
\kwd{BSS}
\kwd{ICA}
\kwd{Jacobi rotations}
\kwd{negentropy}
\kwd{quaternions}
\end{keyword}

% MSC classifications codes, if any
%\begin{keyword}[class=AMS]
%\kwd[Primary ]{}
%\kwd{}
%\kwd[; secondary ]{}
%\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%% start of article main body
% <put your article body there>

%%%%%%%%%%%%%%%%
%% Background %%
%%
\section{Introduction}

Blind source separation (BSS) problems occur in various engineering applications, including the multichannel speech enhancement, image restoration, analysis of electroencephalographic (EEG) signals, and telecommunications \cite{Nai14}.
The BBS consists of recovering of unobservable source signals from their observed mixtures.
Most often, the following linear model \cite{Com94} of the observation vector is assumed:
%
\begin{equation}
{\* x}= {\* A} {\* s},
\label{eq-model}
\end{equation}
%
where $\*A \in \mathbb{R}^{n \times d}$ is an unknown mixing matrix, and $\*s$ is a random vector of $d$ unknown source signals.
Please note that there exist, other models that can deal, for example, with spherical noise \cite{Vir19}, under-determination, and mixture synchronization issues.
In the under-determined case, the number of source signals is greater than the number of mixtures.
Obviously, it makes the BSS problem ill-posed.
In such cases, a concept of sparsity is often employed for a better representation of source signals \cite{Zib01, Li04}.
To tackle the synchronization issues a compressively sensed BSS (CS-BSS) model \cite{Kle12} was proposed which combines a compressive sensing theory \cite{Ran18} with a linear mixing model.
%allows the mixtures to be sampled independently of each other. 

In this work we consider the conventional model \eqref{eq-model} due to its simplicity, and unless otherwise stated, assume that $n = d$.
Thus from a mathematical point of view, we are looking for some matrix $\* W \in \mathbb{R}^{n \times n}$ such that
%
\begin{equation}
{\* y} = {\* W} {\bf x} = {\*D \*P} {\* s},
\end{equation}
%
where $\*D$ and $\*P$ are scaling and permutation matrices, respectively.
This means that the BSS problem can be solved only up to sign, scale, and permutation of the output signals \cite{Hyv00}.
%The presented model is also called

Independent component analysis (ICA) is the most natural approach to solving the BBS problem.
Namely, the ICA aims at transforming multivariate data into components that are as statistically independent from each other as possible.
The ICA can be also viewed as an optimization method that explicitly maximizes some measure of statistical independence between the sources.
Most frequently, it is based on the Central Limit Theorem and attempts to find directions in which some measure of non-Gaussianity is maximized \cite{Hyv99}.
Theoretically, it is possible to extract all components as long as no more than one source is Gaussian and all of the non-Gaussian sources are mutually independent.
In other words, only non-Gaussian components can be separated.

In fact, many ICA algorithms are known.
For general overviews, see for example \cite{Com10, Yu14, Nai14, Cha14}.
The classical ICA methods include the InfoMax \cite{Bel95, Lee99}, FastICA \cite{Hyv99}, and JADE \cite{Car93}.
In this study, we are interested in so-called orthogonal approaches.
Since independence implies uncorrelatedness, these methods directly constrain the unmixed components to be uncorrelated, i.e., for zero-mean and unit-variance signals, it is assumed that $E\{ \*y\*y^T \} = \*I$.
This constraint is usually satisfied by whitening data before rotating them \cite{Hyv99, Car99, Abl18}:
%
\begin{equation}
{\bf y} = {\*B} {\*C}_{\*x\*x}^{-1/2}{\bf x},
\end{equation}
%
where the matrix $\*B$ is constrained to be orthogonal.
The whitening transformation can be computed as follows:
%
\begin{equation}
{\bf C_{xx}}^{-1/2} = {\bf U} {\bf \Lambda}^{-1/2} {\bf U}^T,
\label{eq-wht}
\end{equation}
%
where ${\bf U}$ is the matrix of eigenvectors of ${\bf C_{xx}} = E \{ \*x \*x^T \}$, and ${\bf \Lambda}$ is the diagonal matrix of the corresponding eigenvalues.

Among the orthogonal approaches, the FastICA method \cite{Hyv99} and their variants \cite{Mie19} are probably the most popular techniques.
They are based on maximizing some nonlinear contrast functions that approximate negentropy.
It was shown in \cite{Hyv00} that finding a demixing matrix $\*W$ that minimizes the mutual information \cite{Bel95} is roughly equivalent to finding directions in which the negentropy is maximized.
Therefore, such contrast functions can also be viewed as the measures of non-Gaussianity.
There are two basic variants of the FastICA method: a deflation-based and symmetric algorithm.
In the case of the deflation approach, the independent components are estimated successively one by one and after every iteration step, the Gram-Schmidt orthogonalization is performed.
The major drawback of this approach is that estimation errors of the first vectors are cumulated in the subsequent ones by the orthogonalization \cite{Nor11}.
The symmetric variant of the FastICA algorithm \cite{Hyv99b} estimates the components in parallel.
This consists in computing the one-unit updates for each component, followed by subsequent symmetric orthogonalization of the estimated demixing matrix after each iteration step.

Symmetric orthogonalization in its classical form involves calculating the square root of the inverse matrix in \eqref{eq-wht}, which could be computationally expensive, for large $n$.
Alternatively, simpler iterative methods can be used \cite{Hyv99}.
Nevertheless, the orthogonalization presents some difficulties, especially when the FastICA is implemented in hardware architectures \cite{Pla14}.
Moreover, as we will show, the algorithms with symmetric orthogonalization could be more vulnerable to performance degradation, when the assumption regarding the maximum number of Gaussian sources is not met.

The ICA literature also contains examples of methods that do not require orthogonalization, for instance, Jacobi algorithms \cite{Com94, Car99, Lea03}.
In these approaches, the matrix $\*B$ is considered to be a product of orthogonal matrices, specifically the Jacobi rotations.
The ICA is carried out by optimizing a contrast function for each pair of components that corresponds to a given rotation plane.
Thus, by sweeping the local optimization over all possible pairs, the global independent sources can be extracted.
The most important advantage of the Jacobi algorithms is that the local optimization can be solved analytically if the contrast function is given as a $4$th-order polynomial (e.g. kurtosis).
An example is the JADE technique \cite{Car93} which consists in jointly approximately diagonalizing a set of eigenmatrices.
The eigenmatrices are estimated by $4$th-order moments \cite{Mie15}.
In fact, this is a variant of the Jacobi algorithm, because the diagonalization is performed using the Jacobi rotations.
Unfortunately, the use of the $4$th-order moments as a contrast function may be discouraged due to poor asymptotic efficiency for super-Gaussians and lack of robustness to outliers \cite{Hyv97c}.
Also, the computational complexity of the Jacobi-like algorithms can be relatively high, since the number of required data rotations increases quadratically with $n$.
Theoretically, the number can be reduced by replacing the conventional Jacobi rotations with higher-dimensional rotations.
For instance, $4 \times 4$ orthogonal matrices that represent rotations in $R^4$ can be used.
It has been shown in \cite{Mac95} that such a matrix can be uniquely described by only two unit quaternions.

In our preliminary work \cite{Bor18}, we proposed a novel four-channel ICA algorithm that uses quaternion-based factorizations of $4 \times 4$ orthogonal matrices and exploits a contrast function based on negentropy approximation.
A solution based on the gradient method with the L2-norm constraint was developed.
In this article, we propose to use the Newton-Raphson iterative scheme for solving the optimization task that results in a faster convergence rate as compared to the previous solution.
Moreover, we extended our algorithm for the case of $n \geq 4$ channels, where $n$ is even, by breaking $n$dimensional problem up into a set of four-dimensional subproblems.
For this purpose, a new sweep pattern for data-driven Jacobi-like algorithms is developed.
Also, a more detailed comparative evaluation is performed that involves not only symmetric but also the deflation-based FastICA method and the JADE algorithm.

This paper is organized as follows. Section \ref{sec-quaternions} describes quaternionic factorizations of $4 \times 4$ orthogonal matrices. In Section \ref{sec-method}, a novel ICA algorithm is presented, and the solution based on the Newton-Raphson iterative scheme is derived. Also, a new sweep pattern for data-driven algorithms is developed. Section \ref{sec-experiments} investigates the performance of the proposed method via numerical simulations. Finally the conclusions are given in Section \ref{sec-conclusions}.

\section{Quaternions and rotations in $\mathbb{R}^4$}
\label{sec-quaternions}

Most frequently a quternion $Q \in \mathbb{H}$ is represented in the rectangular form: 	
%
\begin{equation}
Q = q_0 + \textbfit{i}q_1 + \textbfit{j}q_2 + \textbfit{k}q_3, \quad q_0, q_1, q_2, q_3 \in \mathbb{R},
\end{equation}
%
where $\textbfit{i}, \textbfit{j}, \textbfit{k}$ denote imaginary units.
The real part of $Q$ is $q_0$ and the pure quaternion part is $\textbfit{i}q_1 + \textbfit{j}q_2 + \textbfit{k}q_3$.
The multiplication of quaternions is determined by the following rules:
%
\begin{equation}
\textbfit{i}^2 = \textbfit{j}^2 = \textbfit{k}^2 = \textbfit{ijk} = -1.
\end{equation}
%
It is associative and distributes over vector addition, but it is not commutative.
Similarly as for ordinary complex numbers, the conjugate of $Q$ is given by $\bar Q = q_0 - \textbfit{i}q_1 - \textbfit{j}q_2 - \textbfit{k}q_3$, and the norm (modulus), is defined as
%
\begin{equation}
|Q| = \sqrt{q_0^2 + q_1^2 + q_2^2 + q_3^2}.
\end{equation}
%
Since any quaternion $Q$ can be also represented by the vector ${\bf q} = [q_0\, q_1\, q_2\, q_3]^T$, a set $\mathbb{H}$ can be identified with the $\mathbb{R}^4$ vector space. 

Let $\{P, Q\}$ be a pair of unit quaternions i.e. $|P|=|Q|=1$, corresponding to the vectors $\*p, \*q \in \mathbb{R}^4$, and $W$ be a quaternion corresponding to an arbitrary vector $\*w \in \mathbb{R}^4$. 
Consider a real linear transformation from $\mathbb{H}$ to $\mathbb{H}$ that maps $W$ to $PW\bar Q$. 
It can be shown that this transformation can be uniquely represented by a product of the $4 \times 4$ orthogonal matrices \cite{Mac95}:
%
\begin{equation}
{\*M}^{+}(\*p) = 
\begin{bmatrix}[r]
p_0 & -p_1 & -p_2 & -p_3 \\
p_1 &  p_0 & -p_3 &  p_2 \\
p_2 &  p_3 &  p_0 & -p_1 \\
p_3 & -p_2 &  p_1 &  p_0 
\end{bmatrix}, 
\label{eq-M+}
\end{equation}
%
\begin{equation}
{\*M}^{-}(\*q) =
\begin{bmatrix}[r]
 q_0 &  q_1 &  q_2 &  q_3 \\
-q_1 &  q_0 & -q_3 &  q_2 \\
-q_2 &  q_3 &  q_0 & -q_1 \\
-q_3 & -q_2 &  q_1 &  q_0 
\end{bmatrix}.
\label{eq-M-}   
\end{equation} 
%
In particular 
%
\begin{equation}
PW\bar Q \quad \Leftrightarrow \quad {\*M}^{+}(\*p){\*M}^{-}(\*q) \*w. 
\end{equation}
%
It means that for an arbitrary $4 \times 4$ orthogonal matrix $\* R$, there exist a unique pair of unit quaternions $P$ and $Q$ such that $\* R = {\*M}^{+}(\*p){\*M}^{-}(\*q)$.
The proof of this little known theorem can be found in \cite{Bak96}.
It should be noted that this product is commutative, and it is satisfied by the negative quaternions $-P, -Q$ as well.
In fact, the pairs $\{P, Q\}$, and $\{-P, -Q\}$ map to same rotation of $\mathbb{R}^4$.
The possibility of using two parameters to describe a rotation is a special feature of the group of rotations in four dimensions.

\section{Methods}
\label{sec-method}

In this section, we introduce a recently proposed four-channel ICA method \cite{Bor18} and develop a new solution to an optimization problem based on the Newton-Raphson iterative scheme.
We also propose a new ICA algorithm for t$n \geq 4$ channels and give some illustrative examples of possible applications of this approach.

\subsection{Jacobi-based ICA estimation framework} 

Our method strictly follows the Jacobi-based estimation framework.
Namely, the independent components are estimated by recursively applying so-called sweep operations:
%
\begin{equation}
{\*y}_k = {\*R}_k {\*y}_{k-1}, \quad k \geq 1, 
\end{equation}
%
\begin{equation}
{\*y}_0 = {\bf C_{xx}}^{-1/2} {\*x},
\end{equation}
%
where ${\*R}_k$ is an orthogonal matrix and $k$ is a sweep number. 
In the conventional Jacobi-based approaches \cite{Com94, Car99, Lea03}, the matrix ${\*R}_k$ is the product of the rotation matrices:
%
\begin{equation}
   {\bf R}_k = \prod_{i=1}^{n(n-1)/2} {\bf G}(p_i,q_i,\theta_{k,i}),
\label{eq_Rk}
\end{equation}
%
where
%
\begin{equation}
{\bf G}(p,q,\theta) = 
\begin{bmatrix}
 {\bf I}_{p-1} & {\bf 0} & {\bf 0}         & {\bf 0}      & {\bf 0} \\
  {\bf 0} & \cos\theta   & {\bf 0}         & \sin\theta   & {\bf 0} \\
  {\bf 0} & {\bf 0}      & {\bf I}_{q-p-1} & {\bf 0}      & {\bf 0} \\
  {\bf 0} &-\sin\theta   & {\bf 0}         & \cos\theta   & {\bf 0} \\
  {\bf 0} & {\bf 0}      & {\bf 0}         & {\bf 0}      & {\bf I}_{n-q-1} 
\end{bmatrix},
\label{eq-givens}  
\end{equation}
%
represents Jacobi rotation by the $\theta$ angle in the plane determined by the $p$ and $q$ axes. 
Please note that the matrix \eqref{eq-givens} is also known as the Givens matrix and the corresponding Givens rotations are no different from Jacobi rotations \cite{Gol13}.
By using the factorization \eqref{eq_Rk}, an $n$dimensional ICA problem can be decomposed into sequence of $n(n-1)/2$ one-dimensional subproblems. 
Each subproblem consists in searching for the angle that maximizes some contrast function \cite{Car99} in the corresponding rotation plane. 
The data are rotated and the process is repeated cyclically, until all rotation angles converge to zeros.

\subsection{Four-channel algorithm using quaternionic factorization} 

The key idea of the proposed method is to replace the Jacobi rotations with their four-dimensional quaternion-based counterparts.
For simplicty, suppose that $n = 4$. 
Then, six ordinary Jacobi rotations can be replaced with only two quaternion-based rotations as follows:
%
\begin{equation}
\*R_k = \*M^+(\*p_k) \*M^-(\*q_k). 
\label{eq-qfact}
\end{equation}
%
Thus, in each sweep, two vectors $\*p_k$, $\*q_k \in {\mathbb R}^4$ have to be estimated that, for example, maximize some measure of non-Gaussianity of the rotated data. 
We propose to optimize these parameters in two consecutive steps in a similar way as the optimal rotation angles are determined for the Jacobi rotations. 
Firstly, we look for the vector $\*q_k$ and rotate the data using the matrix $\*M^-(\*q_k)$. 
Then we look for $\*p_k$ and rotate the result of the first step by using $\*M^+(\*p_k)$. 
These steps have to be repeated until both rotation matrices are approximately equal to the identity matrices.

It was shown in \cite{Hyv98} that measures of non-Gaussianity based on some non-linear functions are often considerably more robust than cumulant-based approximations.
Moreover, by properly choosing the non-linear function, we can adjust the ICA algorithm to work for practically any distribution of source signals.
One-unit contrast function for measuring non-Gaussianity of any random variable $y$ in the negentropy-based ICA framework is given by \cite{Hyv99}:
%
\begin{equation}
	J_g(y) = [ E \{ g(y) \} - E \{ g(v) \} ] ^2,
\label{eq-Jg}
\end{equation}
%
where $v$ is a standardized Gaussian variable, and $g$ is any non-linear function. 
Examples of such functions can be found in \cite{Tic06}. 
The random variable $y$ is assumed to be zero-mean and unit-variance. 
Please note that the measure \eqref{eq-Jg} is always non-negative, and it equals to zero if $y$ is Gaussian.
As suggested in \cite{Hyv99}, a contrast function for several channels (units) can be obtained by simply maximizing the sum of the one-unit functions.
Thus, by taking into account the unit-norm constraints (that enforces orthogonality), we defined the following optimization problems:
%
\begin{equation}
\begin{gathered}
{\*q}_k = \arg \max_{\*u} \sum_{i=1}^4 
E \{ \bar g([{\*M}^-({\*u})]_{i*} {\*y}^-_k) \}^2 \\
{\rm s.t.} \quad ||{\*u}|| = 1
\label{eq_task1}
\end{gathered}
\end{equation}
%
\begin{equation}
\begin{gathered}
{\*p}_k = \arg \max_{\*u} \sum_{i=1}^4 
E \{ \bar g([{\*M}^+({\*u})]_{i*} {\*y}^+_k) \}^2\\
{\rm s.t.} \quad ||{\*u}|| = 1
\label{eq_task2}
\end{gathered}
\end{equation}
%
where the operator $[{\*A}]_{i*}$ denotes $i$th row of an arbitrary matrix $\bf A$, $\bar g(y) = g(y)- E\{g(v)\}$ is a centered version of the non-linear function $g(y)$, and
%
\begin{equation}
{\*y}^-_k = \*y_{k-1} = \*R_{k-1} \*R_{k-2} \cdots \*R_{1} \*{C_{xx}}^{-1/2} \*x,  \label{eq_rndvec2a}
\end{equation}
\begin{equation}
{\*y}^+_k = {\*M}^-({\*q}_k)  \*y_{k-1}, \label{eq_rndvec2b}
\end{equation}
\begin{equation}
{\*y}_k = {\*M}^+({\*p}_k) {\*M}^-({\*q}_k)  \*y_{k-1},
\end{equation}
% 
are the rotated random vectors for $k \geq 1$.
We assume that the expectations in \eqref{eq_task1}, \eqref{eq_task2} exist and are finite.
Therefore, they can be approximated by sums, and the vector $\*x$ in \eqref{eq_rndvec2a} can be replaced with an observation matrix $\*X$ of size $4 \times m$ consisting of $m$ mixture vectors stacked in column-wise order.
The vectors \eqref{eq_rndvec2a}, \eqref{eq_rndvec2b} can be replaced with transformed data matrices $\*Y_k^-, \*Y_k^+ \in R^{4 \times m}$ as well.
Thus, by using the matrix-vector notation, we can redefine the contrast functions in \eqref{eq_task1} and \eqref{eq_task2} as follows:
%replacing the vectors \eqref{eq_rndvec2a}, \eqref{eq_rndvec2b} with the matrices, say $\*Y^-, \*Y^+ \in R^{4 \times m}$ respectively, and by approximating the expectations in \eqref{eq_task1} and \eqref{eq_task2} by sums, we can redefine maximization problems using the following contrast function:
%
\begin{equation}
J( \*u ) = {\bf 1}_m^T \bar g[ {\*M}({\bf u}) {\*Y_k}]^T \bar g[ {\* M}( \*u ) {\*Y_k}] {\bf 1}_m,
\label{eq-J}
\end{equation}
%
where the terms $\{\*M$, $\*Y_k\}$ stand either for $\{\*M^+$, $\*Y_k^+\}$ or $\{\*M^-$, $\*Y_k^-\}$.
Please note that the problems \eqref{eq_task1} and \eqref{eq_task2} are similar, and the solution depends only on the current sweep data. 
Therefore, in the rest of the article the subscripts $+$, $-$, and the swept number $k$ are dropped to simplify the notation.

\subsection{Solution based on the Newton-Raphson method}

In accordance with the Kuhn-Tucker conditions, any solution of \eqref{eq-J} must satisfy the following equation:
%
\begin{equation}
\nabla_{\*u} J(\*u) + \lambda \*u = \*0,
\label{eq-KKT}
\end{equation}
%
where $\lambda$ is the Lagrange multiplier, and
%
\begin{equation}
\nabla_{\*u} J(\*u) = \left[ 
	\frac{\partial J({\bf u})}{\partial u_0} \cdots
	\frac{\partial J({\bf u})}{\partial u_3} 
\right]^T,
\end{equation}
%
denotes the gradient of $J(\*u)$. 
The partial derivatives can be computed as follows:
%
\begin{equation}
\frac{\partial J({\bf u})}{\partial u_i} = 2 {\bf a}({\bf u})^T {\bf b}_i({\bf u}),
\label{eq-dJdu}
\end{equation}
%
with
\begin{equation}
{\*a}({\*u}) =  \bar g({\bf M}({\bf u}) {\bf Y}) {\bf 1}_m,
\end{equation}
%
\begin{equation}
{\* b}_i({\* u}) = \left[ \bar g'({\bf M}({\*u}) {\*Y}) \circ \frac{\partial {\* M}({\* u})}{\partial u_i} \*Y \right] {\* 1}_m, \quad i = 0,\ldots,3,
\end{equation}
%
where $\circ$ denotes the Hadamard product, $\bar g'$ is the first-order derivative of the non-linear function $\bar g$ and $\frac{\partial {\bf M}({\bf u})}{\partial u_i}$ is $i$th partial derivative of a given quaternionic matrix.
These matrices take the form of simple permutation matrices containing only zeros and ones (see appendix for details), so that, in practice, the multiplication by $\*Y$ can be avoided.
 
One can try to solve the equation \eqref{eq-KKT} numerically by using, for example, the quasi-Newton method in a similar way as in the FastICA approach \cite{Hyv99}.
However, in this case, it could be a more difficult task as it involves a derivation of the second-order derivative (Hessian matrix) of the quadratic form \eqref{eq-dJdu} or finding its approximation.
Even if an explicit expression for this matrix will be obtained, its computation could be very expensive.

In our previous work \cite{Bor18}, we proposed a much simpler solution based on a gradient method  with $L_2$-norm constraint \cite{Dou00}.
A geometric interpretation of this method is given in Fig. \ref{fig-angle}.
In a such approach a pre-scaled gradient of $J(\*u)$, say
%
\begin{equation}
{\*h}(\*u) = \mu \nabla_{\*u} J(\*u),
\label{eq-g} 
\end{equation}
%
with $\mu > 0$ being an empirically chosen scaling factor, is projected onto hyperplane orthogonal to the unit vector $\*u$:
%
\begin{equation}
\boldsymbol \delta = \*h(\*u) - \*u \*u^T \*h(\*u).
\label{eq-delta}
\end{equation}
%
The set of all such vectors at a given point $\*u$ is known as the tangent space of the constraint surface $\| \*u \| = 1$.
In our case, the constraint surface is a unit four-dimensional hypersphere (also known as 3-sphere).
Then, the solution vector is updated by using the following rule: 
%
\begin{equation}
{\*u^*} \leftarrow \*u \cos \theta + \frac{\boldsymbol \delta}{\| \boldsymbol \delta \|} \sin \theta,
\label{eq-update} 
\end{equation}
%
which is equivalent to the rotating $\*u$, the current solution vector, in the direction of $\boldsymbol \delta$ by some angle $\theta$ along the geodesic i.e. curve on the constraint surface that connects two arbitrary points by an arc of shortest length. 

It is a well-known fact that the convergence rate of the gradient methods depends on the value of the scaling factor $\mu$.
In particular, these methods do not necessarily converge for large values of $\mu$.
On the other hand, for small values of $\mu$ the algorithm may converge very slowly.
In the work \cite{Bor18}, we assumed that  $\theta = \| \boldsymbol \delta \|$, and the scaling factor $\mu$ was empirically set to some positive value that depends on a given non-linear function $\bar g(y)$.
Here, we propose a more sophisticated approach that follows the Newton-Raphson iterative scheme and does not require the empirically determined parameters.
Please note that instead of searching for optimal $\mu$, we can equivalently set $\mu = 1$ and maximize the following cost function with respect to rotation angle $\theta$:
%
\begin{equation}
\hat J( \theta ) = {\*1}_m^T  \bar g [ \hat {\*M} (\theta) \*Y]^T \bar g [ \hat {\*M} (\theta) \*Y] {\*1}_m,
\label{eq-Jtheta}
\end{equation}
%
where
%
\begin{equation}
\hat {\bf M}(\theta) = {\bf M} \left( \*u \cos \theta + 
\frac{\boldsymbol \delta}{\|\boldsymbol \delta \|} \sin \theta \right), 
\end{equation}
%
and  $\*u$ is any vector that satisfies $\| {\*u} \| = 1$.
Thus, each constrained optimization task \eqref{eq_task1}, \eqref{eq_task2} can be simplified to one-dimensional angle search on the unit circle.
Now, the Newton-Raphson iterative procedure can be used for finding the optimal rotation angle:
%
\begin{equation}
\theta^* \leftarrow \theta + \left. \frac{\partial \hat J(\theta)}{\partial \theta} \right/ \left| \frac{\partial^2 \hat J(\theta)}{\partial \theta} \right|,
\label{eq-Newton} 
\end{equation}
%
where the first and second derivative of \eqref{eq-Jtheta} are given by
%
\begin{equation}
\frac{\partial \hat J(\theta)}{\partial \theta} = 2 \hat {\*a}(\theta)^T \hat {\*b}(\theta),
\end{equation}
%
\begin{equation}
\frac{\partial^2 \hat J(\theta)}{\partial \theta} = 2 \left[ \hat {\*a}(\theta)^T {\*C}(\theta){\*1}_m,  + \| \hat {\*b}(\theta) \|^2 \right] 
\end{equation}
%
respectively, and
%
\begin{equation}
\hat {\*a}(\theta) = \bar g( \hat {\*M} (\theta) \*Y) {\bf 1}_m,
\label{eq-hat_a}
\end{equation}
%
\begin{equation}
\hat {\*b}(\theta) = \left[ \bar g'( \hat {\*M}(\theta) \*Y) \circ \frac{\partial \hat {\*M}(\theta) }{\partial \theta} \*Y \right] {\*1}_m,
\label{eq-hat_b}
\end{equation}
%
\begin{equation}
{\*C}(\theta) = \bar g'(\hat {\*M}(\theta) \*Y) \circ \left[ \frac{\partial \hat {\*M}(\theta) }{\partial \theta} {\*Y} \right]^{\circ2} - \bar g''(\hat {\*M}(\theta) \*Y) \circ \hat {\*M}(\theta) {\*Y}. 
\label{eq-C}
\end{equation}
%
In the above expression, the operator $\circ2$ denotes element-wise (Hadamard) power, and $\bar g''$ is the second-order derivative of the non-linear function $\bar g$.
Please note that we slightly modified the Newton-Raphson iteration \eqref{eq-Newton} by using plus sign and taking the absolute value of the denominator. 
As suggested in \cite{Mur11} in this way we can enforce the algorithm to find local maxima only. 

In the case of the Jacobi-based algorithms the rotation hyperplanes could overlap each other. 
As a result, local optimizations may not preserve (and often destruct) the global optimality. 
Therefore, during the early sweeps, it may not be necessary to compute the rotation angles with the highest accuracy \cite{Mac95}.
Moreover, each rotation is performed in a local coordinate system.
For these reasons, in each sweep, we perform only one iteration of the Newton-Raphson method per each optimization task \eqref{eq_task1}, \eqref{eq_task2}, assuming that    
%
\begin{equation}
\*u \leftarrow \*u_0 = [1\, 0\, 0\, 0]^T,
\end{equation}
%
\begin{equation}
\theta \leftarrow \theta_0 = 0.
\end{equation}
%
Please note that a similar simplification can also be found in \cite{Oue10}.
It can be easily verified that
\begin{equation}
\frac{\partial \hat {\*M}(\theta) }{\partial \theta} =  {\bf M} \left( \frac{\boldsymbol \delta}{\| \boldsymbol \delta \|} \cos \theta - \*u \sin \theta \right).
\end{equation}
%
Since $\hat {\*M}(0) = \*I$, and $\frac{\partial \hat {\*M}(\theta=0) }{\partial \theta} = {\*M} \left( \frac{\boldsymbol \delta}{\| \boldsymbol \delta \|} \right)$ the expressions \eqref{eq-hat_a} - \eqref{eq-C} can be simplified to:
%
\begin{equation}
\hat {\*a}(0) = \bar g(\*Y){\*1}_m,
\end{equation}
%
\begin{equation}
\hat {\*b}(0) =  \left[ \bar g'(\*Y) \circ \*M \left( \frac{\boldsymbol \delta}{\| \boldsymbol \delta \|} \right) \*Y \right] {\*1}_m,
\end{equation}
\begin{equation}
{\*C}(0) = \bar g'(\*Y) \left[ \*M \left( \frac{\boldsymbol \delta}{\| \boldsymbol \delta \|} \right) \*Y  \right]^{\circ 2} - \bar g''(\*Y) \circ \*Y, 
\end{equation}
%
and the optimal rotation angle can be estimated as follows:
%
\begin{equation}
\theta^* \leftarrow
\frac{\hat {\*a}(0)^T \hat {\*b}(0)}
{|\hat {\*a}(0)^T {\*C}(0){\*1}_m + \| \hat {\*b}(0) \|^2| + \epsilon}.
\label{eq-theta}
\end{equation}
%
where $\epsilon$ is some small positive constant introduced to prevent division by zero.
Once the rotation angle is computed, the current data matrix is transformed as $\*Y \rightarrow \hat {\*M}(\theta^*) \*Y$.
Recalling the factorization \eqref{eq-qfact}, we see that in each sweep, up to two such transformations are needed.
This procedure is repeated cyclically until the convergence is achieved, i.e. for both transformations we have $|\theta^*| < \zeta$, where $\zeta$ is sufficiently small positive constant.

\subsection{Quaternionic ICA for $n \geq 4$ channels}
\label{sec-sweep4d}

An implementation of the four-channel algorithm is rather straightforward as it comprises only two four-dimensional optimization subproblems per sweep.
As shown in the previous section these subproblems can be solved numerically using a quasi-Newton method.
An algorithm for $n \geq 4$ channels can be developed once we determine how to design a sweep pattern consisting of a complete set of four-dimensional subproblems.
Intuitively, every component (row of the data matrix $\*Y$) should be part of any four-dimensional subproblem at least once during the sweep, so no component of the whole data set is exempt from direct participation in the local optimization preserving a drive towards the global optimum.
This idea follows the conventional Jacobi-based ICA estimation framework implemented as a sequence of plane rotations.

It turns out that the sweep patterns comprising four-dimensional subproblems are not difficult to construct and have already been studied in the context of structured eigenproblems \cite{Fab01, Mac05}.
In these approaches, the location of the four-dimensional subproblem was controlled by the position of the element $a_{ij}$ of arbitrary symmetric or skew-symmetric matrix $\*A \in R^{2r \times 2r}$, chosen from the $r \times r$ upper diagonal block, where $r \geq 2$.
This element uniquely determines $4 \times 4$ principal submatrix located in the rows and columns $i, j, r + i$ and $r + j$ as shown below:
%%
%\[\setcellgapes{0pt}
%  \makegapedcells
%\left(
%\begin{array}{ccccc|ccccc}
% 		& \vdots 		& 			& \vdots 		& 			& 		 	& \vdots 			& 			& \vdots 			&  \\
%\cdots 	& a_{i,i}		& \cdots 	& a_{i,j}		& \cdots	& \cdots 	& a_{i,r+i}		& \cdots 	& a_{i,r+j}		& \cdots \\
% 		& \vdots 		& 			& \vdots 		& 			& 		 	& \vdots 			& 			& \vdots 			& \\
%\cdots 	& a_{j,i}		& \cdots 	& a_{j,j}		& \cdots	& \cdots 	& a_{j,r+i}		& \cdots 	& a_{j,r+j}		& \cdots \\
% 		& \vdots 		& 			& \vdots 		& 			& 		 	& \vdots 			& 			& \vdots 			&\\
%\hline 
% 		& \vdots 		& 			& \vdots 		& 			& 		 	& \vdots 			& 			& \vdots 			&  \\
%\cdots 	& a_{r+i,i}	& \cdots 	& a_{r+i,j} 	& \cdots	& \cdots 	& a_{r+i,r+i}		& \cdots 	& a_{r+i,r+j}		& \cdots \\
% 		& \vdots 		& 			& \vdots 		& 			& 		 	& \vdots 			& 			& \vdots 			&  \\
%\cdots 	& a_{r+j,i}	& \cdots 	& a_{r+j,j}	& \cdots	& \cdots 	& a_{r+j,r+i}		& \cdots 	& a_{r+j,r+j}		& \cdots \\
% 		& \vdots 		& 			& \vdots 		& 			& 		 	& \vdots 			& 			& \vdots 			&  \\
%\end{array}\right)
%\]
%
\begin{equation}
\left[ \quad
\begin{matrix}[cccc]
	a_{i,i}	& a_{i,j} & \vrule	& a_{i,r+i} & a_{i,r+j} \\
	a_{j,i}	& a_{j,j} & \vrule	& a_{j,r+i} & a_{j,r+j} \\ \hline
	a_{r+i,i}	& a_{r+i,j} & \vrule	& a_{r+i,r+i} & a_{i,r+j} \\
	a_{r+j,i}	& a_{r+j,j} & \vrule	& a_{r+j,r+i} & a_{j,r+j} \\
\end{matrix}
\quad \right]
\end{equation}
%
where $ 1 \leq i < r$ and $i < j \leq r$.
Please note that such submatrices preserve the structure of the parent matrix $\*A$, i.e. symmetry or skew-symmetry.
Furthermore, a sweep comprised entirely of four-dimensional subproblems can be generated by using any cyclic or quasi-cyclic sweep of the upper-left block of the matrix $\*A$.
As reported in \cite{Mac05}, such a sweep pattern is inherently parallelizable, numerically stable, and shows asymptotic quadratic convergence.

In this work, we propose to extend this idea to a data-driven ICA algorithm for even $n = 2r$.
Instead, to diagonalize the set of eigenmatrices as in the JADE algorithm, we operate directly on the data by following the aforementioned sweep pattern.
Namely, our four-dimensional subproblems are constructed by cyclically selecting the components of indexes $i, j, r + i$, and $r + j$.
All steps of the proposed method are summarized as the pseudocode in Algorithm \ref{alg-QICA}.
In order to generate the pairs $\{i, j\}$ a typical two-dimensional row-cycling ordering was used, but similarly to conventional Jacobi-based methods, other arrangements are also possible \cite{Par20}.

A rigorous theoretical proof of convergence of the algorithm seems to be a challenging task and is out of the scope of this article.
However, we observed empirically that the proposed method is asymptotically quadratically convergent.
Please note that our algorithm requires only $n(n-2)/4$ quaternion-based rotations per sweep in place of $n(n-1)/2$ Jacobi rotations as in the JADE method.
Obviously, it suggests a faster convergence of the proposed method, especially for large $n$.
However, the rotations performed in the JADE method have the lower computational complexity as they are used to transform cumulant matrices only and optimal rotation angles are computed analytically.
These issues will be examined in more detail in the experimental section.
%
%It may be helpful to extend the ideas presented in [*] where the convergence of $4 \times 4$ based quasi-cyclic symmetric Jacobi algorithms was proven.
%As the stopping condition, we propose  
%%
%\begin{equation}
% |\theta^+| + |\theta^-| < \zeta,
%\end{equation}
%%
%where $\zeta$ is a sufficiently small positive scalar constant.
%As reported in \cite{Mac95} in this way it is possible to significantly reduce the number of data rotations and potentially improve the convergence rate for large $n$.
%However, the rotations performed in the JADE method have the low computational complexity as they are used to transform cumulant matrices only and optimal rotation angles are computed analytically.
%These issues will be examinated in more detail in the experimental section.

\begin{algorithm}[tb]
\DontPrintSemicolon
\SetAlgoLined
\KwData{$\*X$ - data matrix of size $n \times m$ containing $n$ observed signals (stacked in rows) each of size $m$ samples, $k_{max}$ - maximum number of sweeps, $\zeta$ - minimum rotation angle}
\KwResult{$\*{\hat Y}$ - matrix of size $n \times m$ containing estimated sources (stacked in rows)}
$\*{C_{xx}} \leftarrow \frac{1}{m} \*X \*X^T$\;
$\*{\hat Y} = [\*{\hat y}_1, \*{\hat y}_2, ..., \*{\hat y}_n]^T \leftarrow \*C_{\*{xx}}^{-1/2} \*X$\;
$\*u_0 \leftarrow [1, 0, 0, 0]^T$ \;
\For {$k \gets 1, 2, ..., k_{max}$}{
	$isRotation \leftarrow false$\;
	\For {$i \gets 1, 2, ..., n/2-1$}{
	 	\For {$j \gets i+1, i+2, ..., n/2$}{
%	 		$\*r \leftarrow \{i, j, n/2+i, n/2+j)\} $ \;
%	 		$\*Y \leftarrow [ \*Y^n]_{\{i,j,n/2+i,n/2+j\}*} $\;
	 		$\*Y \leftarrow [ \*{\hat y}_i, \*{\hat y}_j, \*{\hat y}_{n/2+i}, \*{\hat y}_{n/2+j} ]^T $\;
			compute $\boldsymbol \delta$ using \eqref{eq-delta} with $\*M(\*u) \defeq {\*M}^-(\*u_0), \mu=1$ \;
     			compute $\theta$ angle using \eqref{eq-theta} \;
	 		\If {$|\theta| > \zeta$}{
	 			$\*Y \leftarrow \hat {\*M} (\theta) \*Y$\;
	 			$isRotation \leftarrow true$\;
	 		}
	 		compute $\boldsymbol \delta$ using \eqref{eq-delta} with $\*M(\*u) \defeq {\*M}^+ (\*u_0), \mu=1$ \;
     			compute $\theta$ angle using \eqref{eq-theta} \;
	 		\If {$|\theta| > \zeta$}{
	 			$\*Y \leftarrow \hat {\*M} (\theta) \*Y$\;
	 			$isRotation \leftarrow true$\;
	 		}
	 		$[ \*{\hat y}_i, \*{\hat y}_j, \*{\hat y}_{n/2+i}, \*{\hat y}_{n/2+j} ]^T \leftarrow \*Y$\;
%	 		$[ \*Y^m]_{\{i,j,n/2+i,n/2+j\}*} \leftarrow \*Y$\;
		}
	}	
	\If {$isRotation == false$}{
		break\;
	}
}
\caption{Pseudocode of the proposed method}
\label{alg-QICA}
\end{algorithm}

\subsection{Illustrative examples}

Figure \ref{fig_example1} demonstrates an example of the source recovery using the proposed method.
As input data, a linear mixture of the eight basic signals have been used.
The original and randomly mixed waveforms are depicted in Fig. \ref{fig_example1}a and Fig. \ref{fig_example1}b, respectively.
Figure \ref{fig_example1}c presents the source signals recovered by using the proposed approach.
In order to illustrate the convergence characteristics, in Fig. \ref{fig_example1}d, we also show the negentropy approximation computed after each sweep for all channels as well as for local four-dimensional subproblems.
It can be seen that the method converges in less than six sweeps, and the source signals have been successfully recovered (up to permutation and sign).

We also considered a linear mixture of five images of natural scenes and Gaussian noise.
The test images were selected from the BSDS300 database \cite{Mar01b}.
They were of size $481 \times 321$ pixels with an 8-bit grayscale.
The original images and mixtures are presented in Fig. \ref{fig_example3}a-b.
Figures \ref{fig_example3}c-d present recovered images after $k=2$ and $k=10$ sweeps, respectively.
These experiments show that the proposed algorithm may have good convergence properties in general settings, i.e. for $n > 4$.

In the next example (see Fig. \ref{fig-example2}) we show that the proposed method can also be useful if the number of sources is smaller than four.
Suppose that we have only three independent sources: two speech signals and some noise.
The sources were mixed using a random matrix of size $4 \times 3$ so that there is a linear dependence between the rows of the matrix.
In the case of conventional ICA approaches, the input signals are usually whitened using principal component analysis (PCA), and the dimensionality of the data is decreased at the same time.
Only components that correspond to sufficiently large eigenvalues are taken into account.
However, in the case of our method, we cannot reduce the dimensionality of data as easily as it is closely related to the size of the rotation matrices.
Therefore, for $n<4$ the whitening matrix \eqref{eq-wht} must be properly regularized:
%
\begin{equation}
{\bf C_{xx}}^{-1/2} \approx 
%\begin{bmatrix}[r]
% {\* U}_s & \*0 \\
%\end{bmatrix}
{\* U}
\begin{bmatrix}[c]
 {\* \Lambda}_s^{-1/2} &  \*0 \\
 \*0 &  \*0 \\
\end{bmatrix}
{\* U}^T
%\begin{bmatrix}[r]
% {\* U}_s^T \\
% \*0 \\
%\end{bmatrix}
= {\* U}_s {\* \Lambda_s}^{-1/2} {\* U}_s^T,
\end{equation}
%
where ${\* \Lambda_s}$ is a diagonal matrix of strictly positive eigenvalues of ${\bf C_{xx}}$ and ${\* U}_s$ is a (possibly rectangular) matrix of corresponding eigenvectors.
In other words, the whitening stage is performed in the signal-subspace only, while the minor-components are replaced by zero vectors.
Although such solution is suboptimal from the computational point of view, it makes the implementation much easier.
As can be seen in Fig. \ref{fig-example2}c the original sources have been properly extracted, and the 4th component has been estimated as the zero vector. 
These results suggest that for $n<4$ the minor-components that correspond to nullspace can be treated as non-Gaussian components.
Indeed, the cost function \eqref{eq-Jg} takes positive values for $y$ being a zero signal.
Therefore by maximizing this function we can recover the zero vectors as well as the original sources.

\section{Results and discussion}
\label{sec-experiments}

\subsection{Experimental setup}

The proposed algorithm has been implemented in the MATLAB environment.
For the purpose of empirical evaluation, we denote our method as the QICA (Quaternionic ICA).
It has been compared to the symmetric FastICA, deflation-based FastICA \cite{Hyv99}, and JADE algorithm \cite{Car93}.
For FastICA implementations we used freely available MATLAB packages from \cite{Gav05}.
The implementation of the JADE algorithm is a part of the EEGLAB open-source package \cite{Del04}.
The recent implementation of the JADE algorithm in {R} language as well as some joint diagonalization functions is also described in \cite{Mie17}.

In the case of the QICA and FastICA methods, the experiments were conducted for the three non-linear functions that are frequently considered in publications \cite{Hyv99, Tic06}:
%
\begin{align*}
{\rm POW3:} & \quad \bar g(y) = y^4 - 3,\\
{\rm TANH:} & \quad \bar g(y) = \mathrm{log}( \mathrm{cosh}(y)) - 0.3746,\\
{\rm GAUS:} & \quad \bar g(y) = -\mathrm{exp}(-y^2/2) + 0.7071.
\end{align*}
%
We use these names in order to match the literature, even though one can consider them to be misleading.
Namely, these names are related to the first-order derivatives, not to the functions themselves.

Theoretical stability analysis of the proposed method is out of the scope of this article.
Therefore, the effectiveness of the QICA method was evaluated numerically.
In particular, their convergence characteristics, source separation quality, and computation time were measured for several conditions involving synthetic and real data.
Please note that it makes little sense to measure convergence speed in terms of iteration number or average computation time taken by an algorithm to converge as the reference methods (except the JADE) use different stop criteria.
Therefore, in order to make the comparison meaningful, we decided to disable stop conditions in all methods and perform simulation for a fixed number of iterations measuring the separation quality in each iteration step.
In this way, we were also able to measure the convergence characteristics of the evaluated algorithms.

Synthetic sources were generated for the following distributions: uniform, BPSK, Laplacian, and generalized Gaussians (GG($\alpha$) with parameter $\alpha = 3$, and $0.5$) \cite{Tic06}.
The length of each source signal was $m = 1000$ samples.
In our experiments, we also considered ``speech'' distribution that was represented by data frames randomly selected from real speech recordings.
These recordings were sourced from the publicly available database that was used to develop ITU-T Recommendation \cite{P.50}.
The database consists of short few-seconds sentences spoken by female and male speakers in various languages.
Depending on the scenario, we have selected from 4 up to 12 sentences from the American English set.
Original samples were available in the WAV 32-bit floating-point format with sampling frequency 16kHz, but for our purposes, they were downsampled to 8kHz.

In a basic scenario, we linearly mixed $n = 4$ sources of the same non-Gaussian distribution.
In order to evaluate the robustness of the compared algorithms to the Gaussian sources, we also considered the scenarios in which one or two sources were Gaussian and the others were non-Gaussian but of the same distribution.
Additionally, the experiments with $n=8$ and $n=12$ non-Gaussian sources have been conducted to examine a global convergence of the QICA algorithm.

In all scenarios, the coefficients of the mixing matrix were generated from the uniform distribution. 
Ill-conditioned matrices (with a condition number greater than 1000) were excluded from evaluation.
The performance indexes were averaged over 1000 random realizations of the sources and the mixing matrices.

\subsection{Separation quality}

The separation quality was estimated using signal to interference ratios (SIRs) \cite{Tic06, Tic06b}, and averaged over all non-Gaussian components.
In particular, for the $k$th source signal, the SIR index is given by:
%
\begin{equation}
{\rm SIR}_k = \frac{[\*G]_{kk}^2}{\sum_{l=1}^4[\*G]_{kl}^2},
\end{equation}
%
where $\*G = \hat {\*W} \*A$ is the gain matrix with $\hat {\*W}$ being the estimated demixing matrix of the corresponding algorithm.
Since the BSS problem can be solved only up to permutation, the matrix $\hat {\*W}$ must be appropriately reordered which involves a pairing process \cite{Tic04}.
For this purpose, we followed \cite{Zar10} and used a greedy algorithm.

\begin{table} [b]
\caption{\label{t1} Average SIR (dB) for 4 components with the same distribution, after $k = 40$ iterations.}
\vspace{2mm}
\centerline{
\resizebox{\columnwidth}{!}{
\begin{tabular}{c|ccc|ccc|ccc|c}
 			& \multicolumn{3}{c|}{QICA}	& \multicolumn{3}{c|}{Symm. FastICA}	& \multicolumn{3}{c|}{Deflat. FastICA}	& JADE  \rule{0mm}{3.5mm}\\
PDF		& POW3	& TANH	& GAUSS	& POW3	& TANH	& GAUSS				& POW3	& TANH	& GAUSS			&  \rule{0mm}{3.5mm} \\ \hline
uniform	& 31.49	& 30.44	& 30.36	& 31.50	& 30.45	& 30.37				& 29.36	& 27.90	& 27.80			& \bf31.51 \rule{0mm}{3.5mm} \\
bpsk		& 34.24	& 34.24	& 34.24	& 34.25	& 34.23	& 34.24				& 76.80	& \bf80.88	& 80.33			& 34.24 \rule{0mm}{3.5mm} \\
GG(3.0)	& 21.48	& \bf21.72	& 21.62	& 21.31	& 21.59	& 21.46				& 18.28	& 18.29	& 18.18			& 21.37 \rule{0mm}{3.5mm} \\
Laplacian	& 23.10	& 27.32	& \bf27.57	& 23.48	& 27.40	& \bf27.57				& 20.70	& 24.17	& 24.44			& 23.09 \rule{0mm}{3.5mm} \\
GG(0.5	)	& 27.33	& 31.91	& 31.45	& 27.64	& 31.86	& \bf32.37				& 26.18	& 30.03	& 30.61			& 27.46 \rule{0mm}{3.5mm} \\
speech		& 20.64	& 24.79	& \bf25.58	& 19.65	& 22.82	& 22.98				& 19.13	& 22.06	& 22.67			& 21.68 \rule{0mm}{3.5mm} \\
\end{tabular}}}
\label{tab-1}
\end{table}

\begin{table} [tb]
\caption{\label{t2} Average SIR (dB) for 3 components with the same distribution excluding 1 Gaussian  source, after $k = 40$ iterations.}
\vspace{2mm}
\centerline{
\resizebox{\columnwidth}{!}{
\begin{tabular}{c|ccc|ccc|ccc|c}
 			& \multicolumn{3}{c|}{QICA}	& \multicolumn{3}{c|}{Symm. FastICA}	& \multicolumn{3}{c|}{Deflat. FastICA}	& JADE  \rule{0mm}{3.5mm}\\
PDF		& POW3	& TANH	& GAUSS	& POW3	& TANH	& GAUSS				& POW3	& TANH	& GAUSS			&  \rule{0mm}{3.5mm} \\ \hline
uniform	& \bf29.46	& 28.65	& 28.58	& 26.99	& 27.22	& 27.35				& 28.02	& 26.90	& 26.81			&  29.38 \rule{0mm}{3.5mm} \\
bpsk		& 31.82	& 31.86	& 31.86	& 30.25	& 31.54	& 31.60				& \bf34.18	& 34.12	& 34.12			&  31.83 \rule{0mm}{3.5mm} \\
GG(3.0)	& 19.73	& \bf20.01	& 19.89	& 17.80	& 18.64	& 18.79				& 17.59	& 17.47	& 17.46			&  19.52 \rule{0mm}{3.5mm} \\
Laplacian	& 21.90	& 25.63	& \bf25.98	& 21.92	& 25.05	& 25.30				& 20.26	& 23.62	& 23.98			&  21.89 \rule{0mm}{3.5mm} \\
GG(0.5	)	& 25.21	& 29.61	& \bf30.00	& 25.46	& 29.30	& 29.69				& 24.67	& 28.24	& 28.84			&  25.34 \rule{0mm}{3.5mm} \\
speech		& 15.96	& 21.87	& \bf22.85	& 15.31	& 19.71	& 20.41				& 16.10	& 19.98 	& 20.76			&  18.24 \rule{0mm}{3.5mm} \\
\end{tabular}}}
\label{tab-2}
\end{table}

\begin{table} [tb]
\caption{\label{t3} Average SIR (dB) for 2 components with the same distribution excluding 2 Gaussian  sources, after $k = 40$ iterations.}
\vspace{2mm}
\centerline{
\resizebox{\columnwidth}{!}{
\begin{tabular}{c|ccc|ccc|ccc|c}
 			& \multicolumn{3}{c|}{QICA}	& \multicolumn{3}{c|}{Symm. FastICA}	& \multicolumn{3}{c|}{Deflat. FastICA}	& JADE  \rule{0mm}{3.5mm}\\
PDF		& POW3	& TANH	& GAUSS	& POW3	& TANH	& GAUSS				& POW3	& TANH	& GAUSS			&  \rule{0mm}{3.5mm} \\ \hline
uniform	& 26.89	& 26.34	& 26.23	& 21.74	& 23.84	& 23.93				& 26.50	& 25.58	& 25.44			&  \bf26.97 \rule{0mm}{3.5mm} \\
bpsk		& 29.47	& 29.53	& \bf29.54	& 26.10	& 28.53	& 28.62				& 29.53	& 29.49	& 29.47			&  29.34 \rule{0mm}{3.5mm} \\
GG(3.0)	& 17.54	& \bf18.23	& 18.13	& 14.80	& 16.06	& 16.07				& 16.39	& 16.62	& 16.32			&  17.40 \rule{0mm}{3.5mm} \\
Laplacian	& 20.43	& 23.93	& \bf24.16	& 20.03	& 22.91	& 22.75				& 19.40	& 22.62	& 22.85			&  20.41 \rule{0mm}{3.5mm} \\
GG(0.5	)	& 23.55	& 27.62	& \bf28.03	& 23.72	& 27.06	& 27.21				& 23.28	& 26.87	& 27.31			&  23.69 \rule{0mm}{3.5mm} \\
speech		& 8.85		& 15.75	& \bf16.82	& 7.93		& 13.28	& 13.67				& 8.92		& 14.22	& 15.60			&  9.97 \rule{0mm}{3.5mm} \\
\end{tabular}}}
\label{tab-3}
\end{table}

The results presented in Tables \ref{tab-1}-\ref{tab-3} refer to the scenarios with $n=4$ components containing from 0 up to 2 Gaussian sources.
The maximum number of iterations was set to $k=40$ for all methods.
In the case of the deflation-based FastICA, the components were estimated, one by one, enforcing 40 iterations per component.
The best results obtained for each distribution were denoted using bold fonts.
It can be observed that the performance strongly depends on the distribution of source signals and the number of Gaussian sources.

In the first scenario (Tab. \ref{tab-1} - no Gaussians), we see that the QICA and symmetric FastICA with the POW3 function gives similar results as the JADE algorithm.
It is rather not surprising since the cost function of the JADE algorithm is based on 4th order cumulants and the POW3 function is closely related to the kurtosis.
%Since the SIR indexes were averaged over all components, the symmetric approaches achieved substantially better results than deflation-based FastICA in overall.
%We observed empirically that, for the first component, the performance of the deflation-based FastICA is also similar to that of the other approaches.
%However, for the remaining components, we noticed a substantial performance drop due to deflationary orthogonalization.

In opposition to the JADE algorithm, the performance of the QICA and FastICA methods can be improved for some distributions by selecting appropriate non-linear functions.
It can be observed that if the TANH or GAUSS function is selected, we get higher SIR indexes for mixtures with the Laplacian and GG (0.5) distributions.

By analyzing Tables \ref{tab-1}-\ref{tab-3} it is also clear that the performance of all methods degrades with increasing the number of the Gaussian sources.
However, the QICA method appeared to be more robust in the presence of Gaussian sources.
In the scenario where one Gaussian source was used (Tab. \ref{tab-2}), our algorithm provides substantially better separation quality in almost all cases.
This observation is even more evident in the last scenario (Tab. \ref{tab-3}), where two Gaussians were used.
A similar examination can be made for the JADE method, but the proposed method outperforms the JADE algorithm for the Laplacian and GG(0.5) distributions.

Please note, that if there are no Gaussian sources, then each local minimum or maximum of the one-unit contrast function corresponds to one independent component.
Since the contrast function \eqref{eq-Jg} is blind to Gaussian sources, in the presence of one or more Gaussians, the number of local optima is usually lower than $n$, so that symmetric orthogonalization, as in FastICA approach, may pose some stability issues.
Most frequently, this problem is mitigated by reducing the dimensionality of the observed data, which in fact, involves the identification of non-Gaussian subspace.
The simplest methods for dimensionality reduction assume that information contained in the high-dimensional data is essentially non-Gaussian, and it is located on a lower-dimensional manifold.
Furthermore, the Gaussian components are assumed to be one order of magnitude smaller than the non-Gaussian components.
In such a case, the dimensionality of the observed data can be reduced by using, for example, the PCA technique.
However, this is often a too strict assumption on the data, and in general, the identification of the non-Gaussian subspace is not a trivial task \cite{Bla06, Sas16}.

Unlike the FastICA method, the QICA as well as the JADE approaches are inherently orthogonal and thus more robust to signal model mismatch i.e. situation where there is more than one Gaussian source.
It can be verified empirically that the data rotations are performed in the non-Gaussian subspace only.
Thus, these methods are able to successfully recover non-Gaussian components even if the observation mixtures contain multiple Gaussians.

\subsection{Convergence rate}

Figure \ref{fig-sir} presents average SIR indexes versus iteration number obtained for various data distributions and the number of components ranging from 4 to 12.
In this experiment, no Gaussian sources were used.
In the case of the QICA and symmetric FastICA algorithm, the GAUSS non-linear function was selected.
The deflation-based FastICA is not considered here, as it works differently by recovering the components one by one. 
We observed empirically that, for the first component, the accuracy of this method is similar to that of the other approaches.
However, for the remaining components, we noticed a substantial performance drop due to deflationary orthogonalization.
Thus on average this method performs rather poorly as compared to other algorithms (see Tab. \ref{tab-1}). 

In the case of $n = 4$ channels, we do not observe any significant difference in the convergence rate of the QICA and symmetric FastICA algorithms for most distributions.
However, the proposed method evidently wins with the FastICA, for $n>4$.
The only exception is the BPSK distribution, where the QICA algorithm converges slowly as compared to the reference methods.
Although the SIR decreases with the increase in the number of components, this degradation is similar for all approaches and can be compensated by increasing the number of samples \cite{Tic06}.

It can be observed that the JADE algorithm needs the smallest number of sweeps to converge, but its accuracy is considerably lower than that of the QICA or symmetric FastICA for most distributions.
This is because the JADE algorithm is based on 4th order cumulants which are known to be less sensitive for some mixtures, e.g. Laplacians and super-Gaussians.

Since there is a difference in the number of rotations per sweep between the QICA and JADE algorithms, and both methods can be configured to use the same stop criteria, it may be interesting to measure the total number of rotations that are needed to achieve convergence.
Hence, in the next experiment, the maximum number of sweeps was set to $k_{\max }=100$ and the minimum rotation angle was $\epsilon = 0.01$ for both methods.
Figure \ref{fig-QICAvsJADE} presents the total number of rotations and corresponding SIR index measured for Laplacian distribution versus the number of components.
As indicated in Section \ref{sec-sweep4d}, the QICA method requires approximately 2 times fewer rotations per sweep than the JADE algorithm, but it also needs more sweeps to converge.
Therefore, the total number of rotations is similar for both methods.

\subsection{Computational complexity}

Convergence rate measured by the number of iterations or data rotations, is not the only factor that determines the speed of the algorithms.
In fact, the computational cost of the single iteration, as well as the total execution time, is different for each method.
Figure \ref{fig-times} presents the average execution times measured for various data sizes and algorithms.
Due to significant differences in the implementation of the compared algorithms, we analyze both the execution time of an initialization stage and the time of a single iteration.

It can be seen in Fig. \ref{fig-times}a that the initialization stage of the QICA algorithm is identical to that of the FastICA approach.
Since it only consists in whitening the input data, it is much simpler and less memory demanding than the initialisation stage of the JADE method.
In the case of the JADE method, apart from whitening the input data, it is also necessary to estimate $n(n+1)/2$ cumulant matrices each of size $n \times n$.
Obviously, these operations involve the processing of the whole data set of size $n \times m$ samples, but they are performed only once.
Thus, the memory requirement for the JADE method is of order $O(nm+n^4)$ as opposed to $O(nm)$ for the QICA method.

It can be verified in Fig. \ref{fig-times}b that the execution time of a single sweep of the QICA method increases quadratically with the number of components and linearly with the number of samples. So our method is characterized by the computational complexity of order $O(n^2m)$.
The computation time of a single iteration of the JADE algorithm also increases quadratically with the number of components, but it is constant with respect to the number of samples.
This is because it consists in updating the cumulant matrices only, which does not involve any explicit transformation of the input signals.
Thus its complexity is of order $O(n^2)$

Although the QICA method requires fewer rotations per single sweep than the JADE algorithm (or the same number in total), these rotations are more complex, and similarly to the FastICA method, in every sweep all data samples must be transformed.
As we see in Fig. \ref{fig-times}b, for a relatively small number of samples and a large number of components the average execution time of the single iteration of the proposed method is similar or even smaller than that of the JADE method.
On the other hand, we also see that even for such a small number components as $n=4$ the single iteration time of the QICA method is about six times longer than that of the symmetric FastICA approach.

Although the contrasts for the QICA and FastICA are slightly different, both methods maximize the negentropy measure which involves the evaluation of the same non-linear functions.
Regardless of the implementation, both methods must, at some point, compute these functions and/or their derivatives.
In the case of the QICA approach, for each rotation, three functions: $\bar g$, $\bar g'$, and $\bar g''$ have to be evaluated twice per every four-dimensional subproblem.
Thus at each sweep, we have $3nm(n-2)$ function evaluations in total.
Please note that these operations can be redundant when the rotation angles are approximately close to zero.
The FastICA approach requires evaluating the functions $g'$ and $g''$ only once per component \cite{Hyv99}, so each iteration is related to $2nm$ evaluations in total.

%Therefore, by assuming that $n=4$, single QICA's iteration is at least three times more expensive than the FastICA's iteration.
%As we see in Fig. \ref{fig-times}b this factor may be even greater due to more complicated expressions for the first- and second-order derivatives of the contrasts that are used in the Newton-Raphson procedure. 

The computational complexity of the proposed method can be reduced by taking into account the specific structure of the matrices $\*M^\pm(\*u)$.
In this case, the matrix-vector multiplication can be implemented using only 8 real multiplications \cite{How75}.
On the other hand, multiplication of the $4 \times 4$ matrices can be implemented more efficiently on specific hardware architectures.
Our method consists in decomposing high-dimensional BSS problem into four-dimensional subproblems.
So it well suits modern GPUs and CPUs wchoch support processing four-element vectors.
In particular, typical GPUs are able to read 128 bits in one cycle and one instruction.
Hence, accessing four-element vector-data can be faster than accessing four scalar elements (i.e. 32-bit single-precision floating-point numbers).
Moreover, for some architectures (AMD GPUs), it is essential to organize data into four-element vectors to maximize efficiency \cite{Bua07}.
Please note, that in the case of the CPU-based implementation it is also possible to use 128-bit SSE/AVX instructions.
Other ICA algorithms can be expressed in terms of operations on four-element vectors, but this is unnatural, difficult and results in suboptimal implementations compared to our approach.  

%Although any algorithm can be optimized in this way, in the case of the QICA method such an optimization can be easier to implement as the method inherently breaks large problems up into a set of four-dimensional subproblems.
%In this article, we presented only a four-unit algorithm, but we believe that the QICA method can be extended by breaking the $n$dimensional problem up into a set of four-dimensional subproblems in a similar way as in conventional Jacobi algorithms.
%As reported in \cite{Mac95} in this way it is possible to significantly reduce the number of data rotations and potentially improve the convergence rate for large $n$.
%Please note that our algorithm requires only two four-dimensional rotations per sweep in place of six Jacobi rotations as in the JADE method.
%However, the rotations performed in the JADE method have the low computational complexity as they are used to transform cumulant matrices only and optimal rotation angles are computed analytically.
%Theoretically, it is also possible to modify the QICA method by using cumulant based cost function, but it may result in a performance drop for some distributions.
%Hence, these issues need further investigation.

\section{Conclusion}
\label{sec-conclusions}

We have shown that the quaternionic factorizations of $4 \times 4$ orthogonal matrices can be successfully used in the Jacobi ICA estimation framework for even $n \geq 4$.
Namely, a novel ICA method has been developed, which, on the one hand, follows the Jacobi-based framework, on the other, utilizes the contrast function based on negentropy approximation.
We have shown that the optimization task can be reduced to the problem of a one-dimensional search on an unit circle and solved using the Newton-Raphson method.

The proposed algorithm outperforms the FastICA approach in terms of separation quality if the observation mixtures contain one or more Gaussian sources.
It also shows a faster convergence rate than FastICA approach for most PDFs.
Unlike the JADE algorithm, the proposed method can be adjusted to match a given distribution by selecting an appropriate non-linear function that approximates the negentropy.
The developed method tends to be more computationally demanding than the state of art approaches i.e. symmetric FastICA and JADE algorithm.
However, we believe that the method can be optimized by taking into account specific structure of the quaternion-based rotation matrices.
%using have shown through the experiments that our approach has some advantages over competing methods.

Future works include: rigorous stability analysis of the designed algorithm, further optimizations and developing a GPU-optimized implementation, where the utilization of the four-dimensional rotations seems to be especially promising.

%\appendix
\section*{Appendix }
The derivatives of the quaternionic rotation matrices with respect to the components of the quaternion are given by:

\begin{equation}
\frac{\partial{\*M}^-(\*u)}{\partial u_0} =
%\frac{\partial{\*M}^-}{\partial u_0} =
\begin{bmatrix}[r]
 1 &  0 &  0 &  0 \\
 0 &  1 &  0 &  0 \\
 0 &  0 &  1 &  0 \\
 0 &  0 &  0 &  1 
\end{bmatrix},
\frac{\partial{\*M}^-(\*u)}{\partial u_1} =
%\frac{\partial{\*M}^-}{\partial u_1} =
\begin{bmatrix}[r]
 0 & 1 & 0 & 0 \\
-1 & 0 & 0 & 0 \\
 0 & 0 & 0 &-1 \\
 0 & 0 & 1 & 0 
\end{bmatrix},
\nonumber
\end{equation} 

\begin{equation}
\frac{\partial{\*M}^-(\*u)}{\partial u_2} =
%\frac{\partial{\*M}^-}{\partial u_2} =
\begin{bmatrix}[r]
 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 1 \\
-1 & 0 & 0 & 0 \\
 0 &-1 & 0 & 0 
\end{bmatrix},
\frac{\partial{\*M}^-(\*u)}{\partial u_3} =
%\frac{\partial{\*M}^-}{\partial u_3} =
\begin{bmatrix}[r]
 0 & 0 & 0 & 1 \\
 0 & 0 &-1 & 0 \\
 0 & 1 & 0 & 0 \\
-1 & 0 & 0 & 0 
\end{bmatrix},
\nonumber
\end{equation} 

\begin{equation}
\frac{\partial{\*M}^+(\*u)}{\partial u_0} =
%\frac{\partial{\*M}^+}{\partial u_0} =
\begin{bmatrix}[r]
 1 &  0 &  0 &  0 \\
 0 &  1 &  0 &  0 \\
 0 &  0 &  1 &  0 \\
 0 &  0 &  0 &  1 
\end{bmatrix},
\frac{\partial{\*M}^+(\*u)}{\partial u_1} =
%\frac{\partial{\*M}^+}{\partial u_1} =
\begin{bmatrix}[r]
 0 &-1 & 0 & 0 \\
 1 & 0 & 0 & 0 \\
 0 & 0 & 0 &-1 \\
 0 & 0 & 1 & 0 
\end{bmatrix},
\nonumber
\end{equation} 

\begin{equation}
\frac{\partial{\*M}^+(\*u)}{\partial u_2} =
%\frac{\partial{\*M}^+}{\partial u_2} =
\begin{bmatrix}[r]
 0 & 0 &-1 & 0 \\
 0 & 0 & 0 & 1 \\
 1 & 0 & 0 & 0 \\
 0 &-1 & 0 & 0 
\end{bmatrix},
\frac{\partial{\*M}^+(\*u)}{\partial u_3} =
%\frac{\partial{\*M}^+}{\partial u_3} =
\begin{bmatrix}[r]
 0 & 0 & 0 &-1 \\
 0 & 0 &-1 & 0 \\
 0 & 1 & 0 & 0 \\
 1 & 0 & 0 & 0 
\end{bmatrix}.
\nonumber
\end{equation}  

\section*{List of abbreviations}

\textbf{ICA:}\, Independent component analysis\\
\textbf{JADE:}\, Joint approximate diagonalization of eigenmatrices\\
\textbf{BSS:}\, Blind signal separation\\
\textbf{EEG:}\, Electroencephalography\\
\textbf{PCA:}\, Principal component analysis\\
\textbf{QICA:}\, Quaternionic independent component analysis\\
\textbf{GG:}\, Generalized Gaussian\\
\textbf{SIR:}\, Signal to interference ratio\\
\textbf{GPU:}\, Graphics processing unit\\
\textbf{CPU:}\, Central processing unit\\
\textbf{SSE:}\, Streaming SIMD extensions\\
\textbf{AVX:}\, Advanced vector extensions\\
\textbf{PDF:}\, Probability density function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Backmatter begins here                   %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{backmatter}

\section*{Conset for publication}

Not applicable.

\section*{Availability of data and materials}

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

\section*{Competing interests}

%The author knows of no conflicts of interest associated with this publication.
The authors declare that they have no competing interests.

\section*{Funding}

This work was supported by the Institute of Technical Computer Science and Telecommunications of the Faculty of Computer Science of the Bialystok University of Technology under research subsidy for 2020.

\section*{Authors' contributions}

AB is the sole author of this work and approved the final manuscript.

\section*{Acknowledgements}

Not applicable. 

%\section*{Authors's information (OPTIONAL)}

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if your bibliography is in bibtex format, use those commands:
\bibliographystyle{bmc-mathphys} % Style BST file (bmc-mathphys, vancouver, spbasic).
\bibliography{bmc_article}      % Bibliography file (usually '*.bib' )
% for author-year bibliography (bmc-mathphys or spbasic)
% a) write to bib file (bmc-mathphys only)
% @settings{label, options="nameyear"}
% b) uncomment next line
%\nocite{label}

% or include bibliography directly:
% \begin{thebibliography}
% \bibitem{b1}
% \end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files

\section*{Figures}

% List of figure captions

\begin{figure}[h!]
\centering
%\includegraphics[width=8.5cm]{figure1}
\caption{Optimization using gradient method with $L_2$-norm constraint.}
\label{fig-angle}
\end{figure}

\begin{figure*}[h!]
\centering
%\subfloat[]{\includegraphics[width=8.5cm]{rev_fig2a}}%\hspace{0.5cm}
%\subfloat[]{\includegraphics[width=8.5cm]{rev_fig2b}}\\
%\subfloat[]{\includegraphics[width=8.5cm]{rev_fig2c}}%\hspace{0.5cm}
%\subfloat[]{\includegraphics[width=8.5cm]{rev_fig2d}}
\caption{Demonstration of source extraction using the proposed method. (a) Source signals. (b) Observed mixtures of the source signals. (c) Estimates of the source signals. (d) Contrast function values versus the sweep number.}
\label{fig_example1}
\end{figure*}

\begin{figure*}[h!]
\centering
%\subfloat[]{\includegraphics[width=8.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig3a}}
%\subfloat[]{\includegraphics[width=8.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig3b}}\\
%\subfloat[]{\includegraphics[width=8.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig3c}}
%\subfloat[]{\includegraphics[width=8.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig3d}}
\caption{Separation of images using the proposed method. (a) Original images of natural scenes and Gaussian noise. (b) Linearly mixed images. (c) Separated images after 2 sweeps. (d) Separated images after 10 sweeps.}
\label{fig_example3}
\end{figure*}

\begin{figure*}[h!]
\centering
%\subfloat[]{\includegraphics[width=8.5cm, trim ={0cm 0cm 0cm 0cm}]{figure3a}}\\
%\subfloat[]{\includegraphics[width=8.5cm, trim ={0cm 0cm 0cm 0cm}]{figure3b}}%\hspace{0.5cm}
%\subfloat[]{\includegraphics[width=8.5cm, trim ={0cm 0cm 0cm 0cm}]{figure3c}}
\caption{Extraction of speech signals. (a) Original speech signals. (b) Mixtures of two speech signals and one noise source projected onto four-dimensional space. (c) Source signals estimated using the proposed method.}
\label{fig-example2}
\end{figure*}

\begin{figure*}[h!]
\centering
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig5a}}
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig5b}}
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig5c}}\\
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig5d}}
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig5e}}
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig5f}}
\caption{Average SIR index as a function of sweep/iteration number, obtained using QICA (circles), symmetric FastICA (squares), and JADE algorithm (triangles) for various data distributions. (a) Uniform distribution. (b) BPSK signals. (c) Sub-Gaussian distribution. (d) Laplacian distribution. (e) Super-Gaussian distribution. (f) Speech distribution represented by randomly selected frames of the real speech recordings.}
\label{fig-sir}
\end{figure*}

\begin{figure*}[h!]
\centering
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig6a}}
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig6b}}
\caption{Comparison of the QICA and JADE algorithms for Laplacian distribution and various numbers of components. (a) The number of rotations needed to achieve convergence. (b) SIR index measured after reaching the convergence.}
\label{fig-QICAvsJADE}
\end{figure*}

\begin{figure*}[h!]
\centering
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig7a2}}
%\subfloat[]{\includegraphics[width=5.5cm, trim ={0cm 0cm 0cm 0cm}]{rev_fig7b2}}
\caption{Execution times obtained for QICA, symmetric FastICA, and JADE algorithms for various numbers of samples with fixed $n = 4$ (top) and various numbers of components with fixed $m=1000$ (bottom). (a) Initialization time. (b) Average time of single iteration/sweep.}
\label{fig-times}
\end{figure*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Tables                        %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Use of \listoftables is discouraged.
%%
%\section*{Tables}
%
%\begin{table} [h!]
%\caption{\label{t1} SIR (dB) for 4 components with the same distribution, after $k = 40$ iterations.}
%\vspace{2mm}
%\centerline{
%\resizebox{\columnwidth}{!}{
%\begin{tabular}{c|ccc|ccc|ccc|c}
% 			& \multicolumn{3}{c|}{QICA}	& \multicolumn{3}{c|}{Symm. FastICA}	& \multicolumn{3}{c|}{1-unit FastICA}	& JADE  \rule{0mm}{3.5mm}\\
%PDF		& POW3	& TANH	& GAUSS	& POW3	& TANH	& GAUSS				& POW3	& TANH	& GAUSS			&  \rule{0mm}{3.5mm} \\ \hline
%uniform	& 31.49	& 30.44	& 30.36	& 31.50	& 30.45	& 30.37				& 29.36	& 27.90	& 27.80			& \bf31.51 \rule{0mm}{3.5mm} \\
%bpsk		& 34.24	& 34.24	& 34.24	& 34.25	& 34.23	& 34.24				& 76.80	& \bf80.88	& 80.33			& 34.24 \rule{0mm}{3.5mm} \\
%GG(3.0)	& 21.48	& \bf21.72	& 21.62	& 21.31	& 21.59	& 21.46				& 18.28	& 18.29	& 18.18			& 21.37 \rule{0mm}{3.5mm} \\
%Laplacian	& 23.10	& 27.32	& \bf27.57	& 23.48	& 27.40	& \bf27.57				& 20.70	& 24.17	& 24.44			& 23.09 \rule{0mm}{3.5mm} \\
%GG(0.5	)	& 27.33	& 31.91	& 31.45	& 27.64	& 31.86	& \bf32.37				& 26.18	& 30.03	& 30.61			& 27.46 \rule{0mm}{3.5mm} \\
%\end{tabular}}}
%\label{tab-1}
%\end{table}
%
%\begin{table} [h!]
%\caption{\label{t2} SIR (dB) for 3 components with the same distribution and 1 Gaussian noise source, after $k = 40$ iterations.}
%\vspace{2mm}
%\centerline{
%\resizebox{\columnwidth}{!}{
%\begin{tabular}{c|ccc|ccc|ccc|c}
% 			& \multicolumn{3}{c|}{QICA}	& \multicolumn{3}{c|}{Symm. FastICA}	& \multicolumn{3}{c|}{1-unit FastICA}	& JADE  \rule{0mm}{3.5mm}\\
%PDF		& POW3	& TANH	& GAUSS	& POW3	& TANH	& GAUSS				& POW3	& TANH	& GAUSS			&  \rule{0mm}{3.5mm} \\ \hline
%uniform	& \bf29.46	& 28.65	& 28.58	& 26.99	& 27.22	& 27.35				& 28.02	& 26.90	& 26.81			&  29.38 \rule{0mm}{3.5mm} \\
%bpsk		& 31.82	& 31.86	& 31.86	& 30.25	& 31.54	& 31.60				& \bf34.18	& 34.12	& 34.12			&  31.83 \rule{0mm}{3.5mm} \\
%GG(3.0)	& 19.73	& \bf20.01	& 19.89	& 17.80	& 18.64	& 18.79				& 17.59	& 17.47	& 17.46			&  19.52 \rule{0mm}{3.5mm} \\
%Laplacian	& 21.90	& 25.63	& \bf25.98	& 21.92	& 25.05	& 25.30				& 20.26	& 23.62	& 23.98			&  21.89 \rule{0mm}{3.5mm} \\
%GG(0.5	)	& 25.21	& 29.61	& \bf30.00	& 25.46	& 29.30	& 29.69				& 24.67	& 28.24	& 28.84			&  25.34 \rule{0mm}{3.5mm} \\
%\end{tabular}}}
%\label{tab-2}
%\end{table}
%
%\begin{table} [h!]
%\caption{\label{t3} SIR (dB) for 2 components with the same distribution and 2 Gaussian noise sources, after $k = 40$ iterations.}
%\vspace{2mm}
%\centerline{
%\resizebox{\columnwidth}{!}{
%\begin{tabular}{c|ccc|ccc|ccc|c}
% 			& \multicolumn{3}{c|}{QICA}	& \multicolumn{3}{c|}{Symm. FastICA}	& \multicolumn{3}{c|}{1-unit FastICA}	& JADE  \rule{0mm}{3.5mm}\\
%PDF		& POW3	& TANH	& GAUSS	& POW3	& TANH	& GAUSS				& POW3	& TANH	& GAUSS			&  \rule{0mm}{3.5mm} \\ \hline
%uniform	& 26.89	& 26.34	& 26.23	& 21.74	& 23.84	& 23.93				& 26.50	& 25.58	& 25.44			&  \bf26.97 \rule{0mm}{3.5mm} \\
%bpsk		& 29.47	& 29.53	& \bf29.54	& 26.10	& 28.53	& 28.62				& 29.53	& 29.49	& 29.47			&  29.34 \rule{0mm}{3.5mm} \\
%GG(3.0)	& 17.54	& \bf18.23	& 18.13	& 14.80	& 16.06	& 16.07				& 16.39	& 16.62	& 16.32			&  17.40 \rule{0mm}{3.5mm} \\
%Laplacian	& 20.43	& 23.93	& \bf24.16	& 20.03	& 22.91	& 22.75				& 19.40	& 22.62	& 22.85			&  20.41 \rule{0mm}{3.5mm} \\
%GG(0.5	)	& 23.55	& 27.62	& \bf28.03	& 23.72	& 27.06	& 27.21				& 23.28	& 26.87	& 27.31			&  23.69 \rule{0mm}{3.5mm} \\
%\end{tabular}}}
%\label{tab-3}
%\end{table}

%\begin{table}[h!]
%\caption{Sample table title. This is where the description of the table should go.}
%      \begin{tabular}{cccc}
%        \hline
%           & B1  &B2   & B3\\ \hline
%        A1 & 0.1 & 0.2 & 0.3\\
%        A2 & ... & ..  & .\\
%        A3 & ..  & .   & .\\ \hline
%      \end{tabular}
%\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Additional Files              %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section*{Additional Files}
%  \subsection*{Additional file 1 --- Sample additional file title}
%    Additional file descriptions text (including details of how to
%    view the file, if it is in a non-standard format or the file extension).  This might
%    refer to a multi-page table or a figure.
%
%  \subsection*{Additional file 2 --- Sample additional file title}
%    Additional file descriptions text.


\end{backmatter}
\end{document}
