QR Decomposition of Dual Matrices and its Application to Traveling Wave Identification in the Brain

Renjie Xu Center for Intelligent Multidimensional Data Analysis, Hong Kong Science Park, Shatin, Hong Kong, China. Email: [email protected]. This author is supported by the Hong Kong Innovation and Technology Commission (InnoHK Project CIMDA). Tong Wei Institute of Science and Technology for Brain-Inspired Intelligence, Fudan University, Shanghai, 200433, P.R. China. Email: [email protected]. This author is partially supported by Science and Technology Commission of Shanghai Municipality (No. 23ZR1403000, 20JC1419500, 2018SHZDZX0). Yimin Wei School of Mathematical Sciences and Shanghai Key Laboratory of Contemporary Applied Mathematics, Fudan University, Shanghai, 200433, P.R. China. Email: [email protected]. This author is supported by the National Natural Science Foundation of China under grant 12271108 and the Ministry of Science and Technology of China under grant G2023132005L. Pengpeng Xie Corresponding author (P. Xie). School of Mathematical Sciences, Ocean University of China, Qingdao 266100, P.R. China. [email protected].This author is supported by the National Natural Science Foundation of China under grant 12271108, Fundamental Research Funds for the Central Universities (No. 202264006), Laoshan Laboratory (LSKJ202202302) and Qingdao Natural Science Foundation 23-2-1-158-zyyd-jch. .

Abstract

Matrix decompositions in dual number representations have played an important role in fields such as kinematics and computer graphics in recent years. In this paper, we present a QR decomposition algorithm for dual number matrices, specifically geared towards its application in traveling wave identification, utilizing the concept of proper orthogonal decomposition. When dealing with large-scale problems, we present explicit solutions for the QR, thin QR, and randomized QR decompositions of dual number matrices, along with their respective algorithms with column pivoting. The QR decomposition of dual matrices is an accurate first-order perturbation, with the Q-factor satisfying rigorous perturbation bounds, leading to enhanced orthogonality. In numerical experiments, we discuss the suitability of different QR algorithms when confronted with various large-scale dual matrices, providing their respective domains of applicability. Subsequently, we employed the QR decomposition of dual matrices to compute the DMPGI, thereby attaining results of higher precision. Moreover, we apply the QR decomposition in the context of traveling wave identification, employing the notion of proper orthogonal decomposition to perform a validation analysis of large-scale functional magnetic resonance imaging (fMRI) data for brain functional circuits. Our approach significantly improves the identification of two types of wave signals compared to previous research, providing empirical evidence for cognitive neuroscience theories.

Key words. Dual matrices, QR decomposition, Proper orthogonal decomposition, Randomized algorithm, Traveling wave identification, Brain dynamics

Mathematics Subject Classification. 15A23, 15B33, 65F55

1 Introduction

Time series analysis is a technique for uncovering hidden information and patterns within time-ordered data. Principal Component Analysis (PCA) is employed to interpret the meaning of principal components through orthogonal transformations in time series processes [38]. In statistics, only a small fraction of components contain a significant amount of initial information, and their importance decreases sequentially. However, despite the orthogonality between these principal components, they may not necessarily be independent.

A natural idea is to utilize dual numbers to augment the original data with differentials related to time. This concept has been applied in the context of traveling wave detection [34]. However, this approach exhibits limitations in terms of orthogonality, particularly in the application of brainwave detection. In this paper, we introduce the notion of appropriate orthogonal decomposition, considering the fact that the target time series also represents signal waves. We propose the QR decomposition of dual matrices as an effective and efficient solution to enhance traveling wave detection.

Clifford [5] initially introduced dual numbers in 1873 as a hyper-complex number system, constitute an algebraic structure denoted as $a+b\epsilon$ , where ‘ $a$ ’ and ‘ $b$ ’ are real numbers, and $\epsilon$ satisfies $\epsilon^{2}=0$ with $\epsilon\neq 0$ . These dual numbers can be viewed as an infinitesimal extension of real numbers, representing either minute perturbations of the standard part or time derivatives of time series data.

Dual numbers find applications across diverse domains, including robotics [23], computer graphics [18], control theory [10], optimization [3], automatic differentiation [9], and differential equations [11], owing to their distinctive properties. Consequently, researchers are increasingly exploring the fundamental theoretical aspects of dual matrices, encompassing eigenvalues [21], the singular value decomposition (SVD) [15, 34], the low-rank approximations [19, 36, 37], dual Moore-Penrose pseudo-inverse [1, 6, 17, 27, 28, 29, 31], dual rank decomposition [32], and the QR decomposition [16], among others.

The properties and algorithms for dual matrix decomposition have long been a focal point of research. Previous studies primarily utilized the Gram-Schmidt orthogonalization process or redundant linear systems via Kronecker products to compute matrix decompositions. For instance, Pennestrì et al. [17] employed the Gram-Schmidt orthogonalization process to compute the QR decomposition of dual matrices. Qi et al. [20] used redundant linear systems based on Kronecker products to solve the Jordan decomposition of dual complex matrices with a Jordan block standard part. These approaches are not well-posed for solving large-scale problems.

Recently, Wei et al. [34] introduced a compact dual singular value decomposition (CDSVD) algorithm with explicit solutions. By observing the solving process, we realize that the essence of computing these decompositions is solving a generalized Sylvester equation. Building upon this idea, we provide the QR decomposition of dual matrices (DQR), the thin QR decomposition of dual matrices (tDQR), and the randomized QR decomposition of dual matrices (RDQR), all of which account for column-pivoting scenarios. The QR decomposition of dual matrices also offers an effective approach for computing the dual Moore-Penrose generalized inverse (DMPGI), ensuring heightened precision in the resultant outcomes. The DQR, tDQR, and RDQR proposed here are the proper orthogonal decomposition that we require in time series models. In traveling wave detection, we can see that the DQR, tDQR, and RDQR algorithms perfectly fulfill this function. Moreover, it exhibits stronger orthogonality in large-scale functional magnetic resonance imaging (fMRI) data on brain activity, enabling better identification of traveling wave signals in brain functional areas. In large-scale problems, we can also observe that our proposed QR decomposition is faster than the CDSVD.

The rest of this paper is structured as follows. In Section 2, we present some operational rules and notations for dual numbers and introduce the Sylvester equation. In Section 3, we introduce the concepts of proper orthogonal decomposition for the dual matrix, leading to the development of the QR decomposition (DQR), thin QR decomposition (tDQR), and randomized QR decomposition (RDQR) for the dual matrix, all of which account for column-pivoting scenario. We demonstrate the conditions for the existence of these decompositions and provide algorithms for their computations. In Section 4, we provide a perturbation analysis perspective to demonstrate that the QR decomposition of dual matrices is an accurate first-order perturbation, satisfying the rigorous perturbation bounds for the Q-factor in the QR decomposition. This confirms its improved orthogonality. In Section 5, we compare the computational times of different QR decompositions for large-scale dual matrices and find that their applicability aligns with theoretical expectations. Subsequently, we employed the QR decomposition of dual matrices to compute the DMPGI, thereby attaining results of higher precision. Moreover, we verify that the QR decomposition, based on the idea of proper orthogonal decomposition, can effectively identify traveling waves. Additionally, we demonstrate its utility in accurately characterizing brain functional regions using functional magnetic resonance imaging (fMRI) data on wave patterns.

2 Preliminaries

In this section, we briefly review some basic results for dual matrices and the Sylvester equation. Furthermore, scalars, vectors, and matrices are denoted by lowercase ( $a$ ), bold lowercase ( $\bf{a}$ ), and uppercase letters ( $A$ ), respectively. In some representations of matrix results, we will employ the notation used in MATLAB.

2.1 Dual numbers

In mathematics, dual numbers form an extension of real numbers by introducing a second component with the infinitesimal unit $\epsilon$ satisfying $\epsilon^{2}=0$ and $\epsilon\neq 0$ . This component plays a crucial role in capturing infinitesimal changes and uncertainties in various mathematical operations.

We denote $\mathbb{DR}$ as the set of dual real numbers. Let $d=d_{s}+d_{i}\epsilon$ be a dual real number where $d_{s}\in\mathbb{R}$ and $d_{i}\in\mathbb{R}$ . In this paper, $d_{s}$ and $d_{i}$ are represented as the standard part and the infinitesimal part of $d$ , respectively. If $d_{s}\neq 0$ , we say that $d$ is appreciable. Otherwise, $d$ is infinitesimal. For two real matrices $A_{s},A_{i}\in\mathbb{R}^{m\times n}$ , $A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}$ represents a dual real matrix, where $\mathbb{DR}^{m\times n}$ is the set of dual real matrices. The following proposition shows the preliminary properties of dual real matrices.

Proposition 2.1.

[21] Suppose that $A=(a_{ij})=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}$ , $B=B_{s}+B_{i}\epsilon\in\mathbb{DR}^{n\times p}$ and $C=C_{s}+C_{i}\epsilon\in\mathbb{DR}^{n\times n}$ . Then

1.

The transpose of $A$ is $A^{\top}=(a_{ji})=A_{s}^{\top}+A_{i}^{\top}\epsilon$ . Moreover, $(AB)^{\top}=B^{\top}A^{\top}.$
2.

The dual matrix $C$ is diagonal, if the matrices $C_{s}$ and $C_{i}$ are both diagonal.
3.

The dual matrix $C$ is orthogonal, if $C^{\top}C=I_{n}$ , where $I_{n}$ is an $n\times n$ identity matrix. Furthermore, $CC^{\top}=I_{n}$ , confirms the orthogonal property of the dual matrix $C$ .
4.

The dual matrix $C$ is invertible (nonsingular), if $CD=DC=I_{n}$ for some $D\in\mathbb{DR}^{n\times n}$ . Furthermore, $D$ is unique and denoted by $C^{-1}$ , where $C^{-1}=C_{s}^{-1}-C_{s}^{-1}C_{i}C_{s}^{-1}\epsilon$ .
5.

When $n\geq p$ , the dual matrix $B$ is defined to have orthogonal columns if $B^{\top}B=I_{n}$ . In particular, $B_{s}$ is a real matrix with orthogonal columns, and $B_{s}^{\top}B_{i}+B_{i}^{\top}B_{s}=O_{p}$ , where $O_{p}$ denotes a $p\times p$ zero matrix.

2.2 Sylvester equation

The Sylvester equation, named in honor of the 19th-century mathematician James Joseph Sylvester, represents a prominent matrix equation within the domain of linear algebra. It finds extensive utility across diverse realms of scientific and engineering disciplines, encompassing control theory, signal processing, system stability analysis, and numerical analysis. The Sylvester equation emerges prominently in inquiries linked to time-invariant linear systems, with its solutions serving as indispensable tools for elucidating the characteristics and attributes of these intricate systems. The equation is of the form,

AX+XB=C,

(2.1)

where $A,B$ , and $C$ are given matrices of appropriate dimensions, and $X$ is the unknown matrix that we seek to find. There exists a unique solution $X$ exactly when the spectra of $A$ and $-B$ are disjoint. Here, we first introduce the concept of the Moore-Penrose pseudoinverse.

Definition 2.1.

[13] Given $A\in\mathbb{R}^{m\times n}$ , its Moore-Penrose pseudoinverse, denoted as $A^{\dagger}$ which is a matrix that satisfies the following four properties,

AA^{\dagger}A=A,\quad A^{\dagger}AA^{\dagger}=A^{\dagger},\quad(AA^{\dagger})^{\top}=AA^{\dagger},\quad(A^{\dagger}A)^{\top}=A^{\dagger}A.

(2.2)

Then the solution form of the generalized Sylvester equation can be given by the following theorem.

Theorem 2.1.

[2] Given $A\in\mathbb{R}^{m\times k},B\in\mathbb{R}^{l\times n}$ , and $C\in\mathbb{R}^{m\times n}$ . The equation

AX-YB=C,

(2.3)

has a solution $X\in\mathbb{R}^{k\times n},Y\in\mathbb{R}^{m\times l}$ if and only if

(I_{m}-AA^{\dagger})C(I_{n}-B^{\dagger}B)=O_{m\times n}.

(2.4)

If this is the case, then the general solution of (2.3) has the form

\begin{cases}X=A^{\dagger}C+A^{\dagger}ZB+(I_{k}-A^{\dagger}A)W,\\ Y=-(I_{m}-AA^{\dagger})CB^{\dagger}+Z-(I_{m}-AA^{\dagger})ZBB^{\dagger},\end{cases}

(2.5)

with $W\in\mathbb{R}^{k\times n}$ , $Z\in\mathbb{R}^{m\times l}$ being arbitrary and $A^{\dagger}$ representing the Moore-Penrose pseudoinverse of $A$ .

Theorem 2.1 constitutes a pivotal conclusion, holding significant relevance in the subsequent decomposition of dual matrices.

3 QR Decomposition of Dual Matrices

In prior research [34], the compact SVD decomposition of dual matrices is employed to perform a proper orthogonal decomposition of dual matrices for traveling wave identification. This prompted our exploration of whether QR decomposition could be leveraged for proper orthogonal decomposition in the context of traveling wave identification. In this section, we initially establish the existence of QR decomposition for dual matrices and propose an algorithm. Subsequently, we present an algorithm for QR decomposition with column pivoting in the standard part of dual matrices. We then extend our approach by introducing the concept of thin QR decomposition [13] and presenting algorithms for both thin QR decomposition of dual matrices and thin QR decomposition with column pivoting in the standard part. Furthermore, we harness the idea of randomized QR decomposition with column pivoting [7] to formulate an algorithm for randomized QR decomposition with column pivoting in the standard part of dual matrices. The latter two categories of algorithms are specifically designed for addressing large-scale data scenarios where $m\gg n$ , enables a proper orthogonal decomposition.

3.1 QR Decomposition of Dual Matrices (DQR)

Theorem 3.1.

(DQR Decomposition) Given $A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}$ . Assume that $A_{s}=Q_{s}R_{s}$ is the QR decomposition of $A_{s}$ , where $Q_{s}\in\mathbb{R}^{m\times m}$ is an orthogonal matrix and $R\in\mathbb{R}^{m\times n}$ is an upper triangular matrix. Then the QR decomposition of the dual matrix $A=QR$ where $Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times m}$ is an orthogonal dual matrix and $R=R_{s}+R_{i}\epsilon\in\mathbb{R}^{m\times n}$ is an upper triangular dual matrix exists. Additionally, the expressions of $Q_{i}$ , $R_{i}$ , the infinitesimal parts of $Q,R$ , are

\begin{cases}Q_{i}&=Q_{s}P\\ R_{i}&=Q_{s}^{\top}A_{i}-PR_{s},\end{cases}

(3.1)

where $P\in\mathbb{R}^{m\times m}$ is a skew-symmetric matrix.

Assuming $\operatorname{rank}(A_{s})=K\leq\min(m,n)$ and $R_{s}=(r_{s_{ij}})\in\mathbb{R}^{m\times n}$ , we can derive the lower triangular portion’s expression for $P$ as

\begin{cases}p_{1}(2:m)&=b_{1}(2:m)/r_{s_{11}},\\ p_{2}(3:m)&=(b_{2}(3:m)-r_{s_{12}}p_{1}(3:m))/r_{s_{22}},\\ &\ldots\\ p_{K}(K+1:m)&=(b_{K}(K+1:m)-\sum\limits_{t=1}^{K-1}r_{s_{tK}}p_{t}(K+1:m))/r_{s_{KK}},\end{cases}

(3.2)

where $Q_{s}^{\top}A_{i}=B=[b_{1},b_{2},\ldots,b_{n}]$ and the remaining lower triangular portion is entirely composed of zeros.

Proof.

Assume the QR decomposition of $A$ exists, then we have that

\begin{cases}A_{s}=Q_{s}R_{s}\\ A_{i}=Q_{s}R_{i}+Q_{i}R_{s}\end{cases}

(3.3)

by simple computation.

By treating $R_{i}$ as $X$ and $Q_{i}$ as $Y$ , the given equation $A_{i}=Q_{s}R_{i}+Q_{i}R_{s}$ in (3.3) transforms into a Sylvester equation involving $Q_{i}$ and $R_{i}$ . Leveraging Theorem 2.1, if and only if

(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}(I_{m}-R_{s}^{\dagger}R_{s})=O_{m\times n},

(3.4)

which holds true identically by the orthogonal nature of $Q_{s}$ and the fact $Q_{s}^{\dagger}=Q_{s}^{\top}$ , one can derive the solution

\begin{cases}R_{i}&=Q_{s}^{\dagger}A_{i}+Q_{s}^{\dagger}Z(-R_{s})+(I_{m}-Q_{s}^{\dagger}Q_{s})W\\ Q_{i}&=-(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}(-R_{s})^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\dagger})Z(-R_{s})(-R_{s})^{\dagger}\end{cases}

(3.5)

where $Z\in\mathbb{R}^{m\times m}$ and $W\in\mathbb{R}^{m\times n}$ are arbitrary.

Upon examining the first expression of $R_{i}$ in (3.5), one can deduce, based on the established fact that $Q_{s}^{\dagger}=Q_{s}^{\top}$ ,

$\displaystyle R_{i}$	$\displaystyle=Q_{s}^{\dagger}A_{i}+Q_{s}^{\dagger}Z(-R_{s})+(I_{m}-Q_{s}^{\dagger}Q_{s})W$	(3.6)
	$\displaystyle=Q_{s}^{\top}A_{i}+Q_{s}^{\top}Z(-R_{s})+(I_{m}-Q_{s}^{\top}Q_{s})W$
	$\displaystyle=Q_{s}^{\top}A_{i}-Q_{s}^{\top}ZR_{s}.$

Likewise, drawing on the foundational premise of $Q_{s}^{\dagger}=Q_{s}^{\top}$ and $I_{m}-Q_{s}Q_{s}^{\top}=O_{m}$ , we can derive the second expression of $Q_{i}$ in (3.5) into

$\displaystyle Q_{i}$	$\displaystyle=-(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}(-R_{s})^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\dagger})Z(-R_{s})(-R_{s})^{\dagger}$	(3.7)
	$\displaystyle=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\top})ZR_{s}R_{s}^{\dagger}$
	$\displaystyle=Z.$

Subsequently, the determination of arbitrary values for $Z$ becomes imperative. Notably, the orthogonal nature of the dual matrix $Q$ and the upper triangular configuration of the dual matrix $R$ can be employed as constraining equations to ascertain the values of $Z$ . Primarily, our scrutiny revolves around the orthogonal property of the dual matrix $Q$ . This leads us to consider the following

Q_{i}^{\top}Q_{s}+Q_{s}^{\top}Q_{i}=O_{m}.

(3.8)

Substituting (3.7) into the term $Q_{s}^{\top}Q_{i}$ in (3.8) yields

Q_{s}^{\top}Q_{i}=Q_{s}^{\top}Z.

(3.9)

Hence, (3.8) can be transformed into

Z^{\top}Q_{s}+Q_{s}^{\top}Z=O_{m}.

(3.10)

Here, owing to the orthogonal and invertible nature of matrix $Q_{s}$ , we can denote $P=Q_{s}^{\top}Z\in\mathbb{R}^{m\times m}$ (thus equivalently, $Z=Q_{s}P$ , then the investigation into the values of $Z$ can be transmuted into an exploration of the values pertaining to $P$ ). According to (3.9), it follows that $P$ is a skew-symmetric matrix. Through the skew-symmetric matrix $P=Q_{s}^{\top}Z$ , we can attain expressions for $Q_{i}$ and $R_{i}$ ,

\begin{cases}Q_{i}&=Q_{s}P\\ R_{i}&=Q_{s}^{\top}A_{i}-PR_{s}.\end{cases}

(3.11)

Now, we must revisit the upper triangular property of the dual matrix $R$ to determine $P$ . Upon substituting the expression $P=Q_{s}^{\top}Z$ into (3.6), we obtain

PR_{s}=Q_{s}^{\top}A_{i}-R_{i}.

(3.12)

We denote $B=Q_{s}^{\top}A_{i}\in\mathbb{R}^{m\times n}$ and $B=[b_{1},b_{2},\ldots,b_{n}],R_{i}=[r_{i1},r_{i2},\ldots,r_{in}],R_{s}=[r_{s1},r_{s2},\ldots,r_{sn}]$ where $b_{k},r_{ik},r_{sk}\in\mathbb{R}^{m},k=1,2,\ldots,n$ . By partitioning the matrices column-wise, (3.12) is transformed into the following system of linear equations,

\begin{cases}Pr_{s1}&=b_{1}-r_{i1},\\ Pr_{s2}&=b_{2}-r_{i2},\\ &\ldots\\ Pr_{sk}&=b_{k}-r_{ik},\\ &\ldots\end{cases}

(3.13)

Exploiting the upper triangular nature of matrix $R_{s}$ , along with $R_{s}=(r_{s_{ij}})\in\mathbb{R}^{m\times n}$ and $P=[p_{1},p_{2},\ldots,p_{n}]$ where $p_{k}\in\mathbb{R}^{m},k=1,2,\ldots,m$ , the system (3.13) can be transmuted into

\begin{cases}r_{s_{11}}p_{1}&=b_{1}-r_{i1},\\ r_{s_{12}}p_{1}+r_{s_{22}}p_{2}&=b_{2}-r_{i2},\\ &\ldots\\ r_{s_{1k}}p_{1}+r_{s_{2k}}p_{2}+\ldots+r_{s_{kk}}p_{k}&=b_{k}-r_{ik},\\ &\ldots\end{cases}

(3.14)

Recognizing the upper triangular structure of $R_{i}$ , it follows that for $r_{ik},k=1,2,\ldots,n$ , the lower $n-k$ entries $r_{ik}(k+1:m)$ are all zeros. Consequently, $p_{k}$ also requires computation solely for the lower $n-k$ entries $p_{k}(k+1:m)$ , thereby ensuring the maintenance of $P$ ’s skew-symmetric property.

Given $\operatorname{rank}(R_{s})=\operatorname{rank}(A_{s})=K\leq\min(m,n)$ , it follows that all the elements $r_{s_{ij}}=0,i>K$ . Hence, the system of $K$ equations has the solutions

\begin{cases}p_{1}(2:m)&=b_{1}(2:m)/r_{s_{11}},\\ p_{2}(3:m)&=(b_{2}(3:m)-r_{s_{12}}p_{1}(3:m))/r_{s_{22}},\\ &\ldots\\ p_{K}(K+1:m)&=(b_{K}(K+1:m)-\sum\limits_{t=1}^{K-1}r_{s_{tK}}p_{t}(K+1:m))/r_{s_{KK}}.\end{cases}

(3.15)

Thus, we complete the proof.

Theorem 3.1 furnishes the proof of existence for the QR decomposition of the dual matrices, accompanied by a delineation of the construction process within the proof. In Algorithm 1, we furnish pseudocode for the QR decomposition of the dual matrices.

Input:

A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}

2Compute the QR decomposition of

A_{s}=Q_{s}R_{s}(Q_{s}\in\mathbb{R}^{m\times m},R_{s}\in\mathbb{R}^{m\times n})

3Compute

K=\operatorname{rank}(R_{s})

4Set

P=zeros(m,m)

and

B=Q_{s}^{\top}A_{i}

5Set

i=1

P(2:m,1)=B(2:m,1)/R_{s}(1,1)

7while $i<K$ do

i=i+1

P(i+1:m,i)=(B(i+1:m,i)-\sum\limits_{t=1}^{i-1}R_{s}(t,i)P(i+1:m,t)/R_{s}(i,i).

P(i,i+1:m)=-P(i+1:m,i)^{\top}

13Set

Q_{i}=Q_{s}P

14Set

R_{i}=B-PR_{s}

Output:

Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times m}

is an orthogonal dual matrix, and

R=R_{s}+R_{i}\epsilon\in\mathbb{DR}^{m\times n}

is an upper triangular dual matrix.

Algorithm 1 QR Decomposition of Dual Matrices (DQR)

In practical traveling wave detection [34], it becomes necessary to perform column pivoting in the standard components of the QR decomposition. Thus, suppose that the QR decomposition with column pivoting in the standard part of $A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}$ exists, then we have

\begin{cases}A_{s}P_{s}&=Q_{s}R_{s}\\ A_{i}P_{s}&=Q_{s}R_{i}+Q_{i}R_{s}.\end{cases}

(3.16)

Building upon (3.16), we can employ a similar approach as in Theorem 3.1 to complete the proof. We omit the details here. Algorithm 2 presents the pseudocode for the QR decomposition with column pivoting in the standard part of the dual matrix (DQRCP).

Input:

A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}

2Compute the QR decomposition with column pivoting of

A_{s}P_{s}=Q_{s}R_{s}(Q_{s}\in\mathbb{R}^{m\times m},R_{s}\in\mathbb{R}^{m\times n},P_{s}\in\mathbb{R}^{n\times n})

3Set

P=zeros(m,m)

and

B=Q_{s}^{\top}A_{i}P_{s}

4Set

i=1

P(2:m,1)=B(2:m,1)/R_{s}(1,1)

6while $i<K$ do

i=i+1

P(i+1:m,i)=(B(i+1:m,i)-\sum\limits_{t=1}^{i-1}R_{s}(t,i)P(i+1:m,t)/R_{s}(i,i)

P(i,i+1:m)=-P(i+1:m,i)^{\top}

12Set

Q_{i}=Q_{s}P

13Set

R_{i}=B-PR_{s}

Output:

Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times m}

is an orthogonal dual matrix,

R=R_{s}+R_{i}\epsilon\in\mathbb{DR}^{m\times n}

is an upper triangular dual matrix, and

P_{s}\in\mathbb{R}^{n\times n}

is a permutation matrix.

Algorithm 2 QR Decomposition with Column Pivoting in the standard part of Dual Matrices (DQRCP)

3.2 Thin QR Decomposition of Dual Matrices (tDQR)

In Algorithm 1, it is evident that while computing the QR decomposition of the dual matrix $A\in\mathbb{R}^{m\times n}$ , in order to construct the skew-symmetric matrix $P$ , a space of dimension $m^{2}$ is required. In the context of large-scale numerical computations, a large value of $m$ can have a detrimental impact on our calculations. Therefore, this necessitates our contemplation of thin QR decomposition of dual matrices. The concept of thin QR decomposition was initially introduced by Golub and Van Loan [13]. For a matrix $A\in\mathbb{R}^{m\times n}(m\geq n)$ with full column rank, the thin QR decomposition $A=QR$ entails the unique existence of a matrix $Q\in\mathbb{R}^{m\times n}$ with columns orthogonality and an upper triangular matrix $R\in\mathbb{R}^{n\times n}$ with positive diagonal entries. Subsequently, we proceed to present the thin QR decomposition of the dual matrices (tDQR decomposition) based on this notion.

Theorem 3.2.

(tDQR Decomposition) Given $A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}(m\geq n)$ where $A_{s}\in\mathbb{R}^{m\times n}$ is full column rank. Assume that $A_{s}=Q_{s}R_{s}$ is the thin QR decomposition of $A_{s}$ , where $Q_{s}\in\mathbb{R}^{m\times n}$ is a matrix with orthogonal columns and $R\in\mathbb{R}^{n\times n}$ is an upper triangular matrix with positive diagonal entries. Then the thin QR decomposition of the dual matrix $A=QR$ where $Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times n}$ is a dual matrix with orthogonal columns and $R=R_{s}+R_{i}\epsilon\in\mathbb{R}^{n\times n}$ is an upper triangular dual matrix exists. Additionally, the expressions of $Q_{i}$ , $R_{i}$ , the infinitesimal parts of $Q,R$ , are

\begin{cases}Q_{i}&=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}+Q_{s}P\\ R_{i}&=Q_{s}^{\top}A_{i}-PR_{s}\end{cases}

(3.17)

where $P\in\mathbb{R}^{n\times n}$ is a skew-symmetric matrix.

We can derive the lower triangular portion’s expression for $P$ as

\begin{cases}p_{1}(2:n)&=b_{1}(2:n)/r_{s_{11}},\\ p_{2}(3:n)&=(b_{2}(3:n)-r_{s_{12}}p_{1}(3:n))/r_{s_{22}},\\ &\ldots\\ p_{n-1}(n)&=(b_{n-1}(n)-\sum\limits_{t=1}^{n-2}r_{s_{t(n-1)}}p_{t}(n))/r_{s_{(n-1)(n-1)}},\end{cases}

(3.18)

where $Q_{s}^{\top}A_{i}=B=[b_{1},b_{2},\ldots,b_{n}]$ .

The proof can be found in the Appendix A.

Theorem 3.2 provides the thin QR decomposition of dual matrices (tDQR), offering an efficient alternative for the QR decomposition of dual matrices (DQR) when dealing with large-scale problems where $m\gg n$ . Algorithm 3 furnishes pseudocode for the tDQR.

Input:

A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}

where

A_{s}

is of full column rank.

2Compute the thin QR decomposition of

A_{s}=Q_{s}R_{s}(Q_{s}\in\mathbb{R}^{m\times n},R_{s}\in\mathbb{R}^{n\times n})

3Set

P=zeros(n,n)

and

B=Q_{s}^{\top}A_{i}

4Set

i=1

P(2:n,1)=B(2:n,1)/R_{s}(1,1)

6while $i<n-1$ do

i=i+1

P(i+1:n,i)=(B(i+1:n,i)-\sum\limits_{t=1}^{i-1}R_{s}(t,i)P(i+1:n,t)/R_{s}(i,i)

P(i,i+1:n)=-P(i+1:n,i)^{\top}

12Set

Q_{i}=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}+Q_{s}P

13Set

R_{i}=B-PR_{s}

Output:

Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times n}

is a dual matrix with orthogonal columns, and

R=R_{s}+R_{i}\epsilon\in\mathbb{DR}^{n\times n}

is an upper triangular dual matrix.

Algorithm 3 Thin QR Decomposition of Dual Matrices (tDQR)

Similar to Algorithm 2, we can derive an algorithm for the thin QR decomposition with column pivoting in the standard part of the dual matrix based on analogous reasoning. Algorithm 4 presents the pseudocode for the thin QR decomposition with column pivoting in the standard part of the dual matrix (tDQRCP).

Input:

A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}

, where

A_{s}

is of full column rank.

2Compute the thin QR decomposition with column pivoting of

A_{s}P_{s}=Q_{s}R_{s}(Q_{s}\in\mathbb{R}^{m\times n},R_{s}\in\mathbb{R}^{n\times n},P_{s}\in\mathbb{R}^{n\times n})

3Set

P=zeros(n,n)

and

B=Q_{s}^{\top}A_{i}P_{s}

4Set

i=1

P(2:n,1)=B(2:n,1)/R_{s}(1,1)

6while $i<n-1$ do

i=i+1

P(i+1:n,i)=(B(i+1:n,i)-\sum\limits_{t=1}^{i-1}R_{s}(t,i)P(i+1:n,t)/R_{s}(i,i)

P(i,i+1:n)=-P(i+1:n,i)^{\top}

12Set

Q_{i}=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}R_{s}^{-1}+Q_{s}P

13Set

R_{i}=B-PR_{s}

Output:

Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times m}

is a dual matrix with orthogonal columns,

R=R_{s}+R_{i}\epsilon\in\mathbb{DR}^{m\times n}

is an upper triangular dual matrix, and

P_{s}\in\mathbb{R}^{n\times n}

is a permutation matrix.

Algorithm 4 Thin QR Decomposition with Column Pivoting in the standard part of Dual Matrices (tDQRCP)

3.3 Randomized QR Decomposition with Column Pivoting of Dual Matrices (RDQRCP)

In practical applications utilizing QR decomposition with column pivoting, it is often sufficient to compute only the first $k$ rows of $R$ . Given an $m\times n$ matrix $A$ , we seek a truncated QR decomposition for $k\leq\min(m,n)$

AP\approx QR,

(3.19)

where $Q\in\mathbb{R}^{m\times k}$ is orthonormal, $P\in\mathbb{R}^{n\times n}$ is a permutation, and $R\in\mathbb{R}^{k\times n}$ is upper triangular. Addressing this concern, Duersch and Gu devised a randomized algorithm for QR decomposition with column pivoting [7]. This truncated QR decomposition employing randomized algorithms proves to be highly effective for numerous low-rank matrix approximation techniques. According to the results of the Johnson-Lindenstrauss lemma [33], randomized sampling algorithms maintain a high probability guarantee of approximation error while simultaneously reducing communication complexity through dimensionality reduction. The Johnson-Lindenstrauss lemma implies that let $a_{j}$ represent the $j$ -th column ( $j=1,2,\ldots,n$ ) of $A$ . We construct the sample columns $b_{j}=\Omega a_{j}$ by a Gaussian Independent Identically Distributed (GIID) compression matrix $\Omega\in\mathbb{R}^{l\times m}$ . The 2-norm expectation values and variance of $b_{j}$ are,

\mathbb{E}(\|b_{j}\|_{2}^{2})=l\|a_{j}\|_{2}^{2}\text{ and }\mathbb{V}(\|b_{j}\|_{2}^{2})=2l\|a_{j}\|_{2}^{4}

(3.20)

Moreover, the probability of successfully detecting all column norms as well as all distances between columns within a relative error $\tau$ for $0<\tau<\frac{1}{2}$ is bounded by

P\left(\left|\frac{\|b_{j}-b_{i}\|_{2}^{2}}{l\|a_{j}-a_{i}\|_{2}^{2}}-1\right|\leq\tau\right)\geq 1-2e^{\frac{-l\tau^{2}}{4}(1-\tau)}.

(3.21)

We can select the column with the largest approximate norm by using the sample matrix $B=\Omega A$ . The remaining columns can be processed as represented in Algorithm 5, with detailed proofs available in the work by Duersch and Gu [7].

Input:

A\in\mathbb{R}^{m\times n}

k\ll\min(m,n)

the desired approximation rank.

2Set sample rank

l=k+p

needed for acceptable sample error.

3Generate random

l\times m

GIID compression matrix

\Omega

and form the sample

B=\Omega A

4Get the QR decomposition with column pivoting of the sample,

BP=Q_{b}R_{b}(Q_{b}\in\mathbb{R}^{l\times l},R_{b}\in\mathbb{R}^{l\times n},P\in\mathbb{R}^{n\times n})

5Apply permutation

A^{(1)}=AP

6Compute the thin QR decomposition of

A^{(1)}(:,1:k)=QR_{11}(Q\in\mathbb{R}^{m\times k},R_{11}\in\mathbb{R}^{k\times k})

7Finish

k

rows of

R

in remaining columns,

R_{12}=Q^{\top}(:,1:k)A^{(1)}(:,k+1:end)

Output:

Q\in\mathbb{R}^{m\times k}

is a matrix with orthogonal columns,

R\in\mathbb{R}^{k\times n}

is a truncated upper trapezoidal matrix, and

P\in\mathbb{R}^{n\times n}

is a permutation matrix.

Algorithm 5 Single-Sample Randomized QRCP (RQRCP)

Building upon the randomized QR decomposition with column pivoting for (3.19) as illustrated in Algorithm 5, we present a randomized QR decomposition with column pivoting in the standard part of dual matrices.

Theorem 3.3.

(RDQRCP Decomposition) Given $A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}(m\geq n)$ . Assume that $A_{s}P_{s}=Q_{s}R_{s}$ is the randomized QR decomposition with column pivoting of $A_{s}$ , where $k\leq\min(m,n)$ , $Q_{s}\in\mathbb{R}^{m\times k}$ is a matrix with orthogonal columns, $R\in\mathbb{R}^{k\times n}$ is a truncated upper trapezoidal matrix and $P_{s}\in\mathbb{R}^{n\times n}$ is a permutation matrix and $\operatorname{rank}(A_{s})\geq k$ . Then the randomized QR decomposition with column pivoting in the standard part of the dual matrix $AP_{s}=QR$ where $Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times k}$ is a dual matrix with orthogonal columns and $R=R_{s}+R_{i}\epsilon\in\mathbb{R}^{k\times n}$ is a truncated upper trapezoidal dual matrix exists if and only if

(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}(I_{n}-R_{s}^{\dagger}R_{s})=O_{m\times n}.

(3.22)

Additionally, the expressions of $Q_{i}$ , $R_{i}$ , the infinitesimal parts of $Q,R$ , are

\begin{cases}Q_{i}&=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}R_{s}^{\dagger}+Q_{s}P\\ R_{i}&=Q_{s}^{\top}A_{i}P_{s}-PR_{s},\end{cases}

(3.23)

where $P\in\mathbb{R}^{k\times k}$ is a skew-symmetric matrix.

We can derive the lower triangular portion’s expression for $P$ as

\begin{cases}p_{1}(2:k)&=b_{1}(2:k)/r_{s_{11}},\\ p_{2}(3:k)&=(b_{2}(3:k)-r_{s_{12}}p_{1}(3:k))/r_{s_{22}},\\ &\ldots\\ p_{k-1}(k)&=(b_{k-1}(k)-\sum\limits_{t=1}^{k-2}r_{s_{t(k-1)}}p_{t}(k))/r_{s_{(k-1)(k-1)}},\end{cases}

(3.24)

where $Q_{s}^{\top}A_{i}P_{s}=B=[b_{1},b_{2},\ldots,b_{k},\ldots,b_{n}]$ .

The proof can be found in the Appendix A.

Theorem 3.3 provides the existence of a randomized QR decomposition with column pivoting in the standard part of dual matrices. Algorithm 6 presents its pseudocode.

Input:

A=A_{s}+A_{i}\epsilon\in\mathbb{DR}^{m\times n}

2Compute the RQRCP of

A_{s}P_{s}=Q_{s}R_{s}

in Algorithm 5

(Q_{s}\in\mathbb{R}^{m\times k},R_{s}\in\mathbb{R}^{k\times n},P_{s}\in\mathbb{R}^{n\times n})

3while $(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}(I_{n}-R_{s}^{\dagger}R_{s})=O_{m\times n}$ do

5 Set

P=zeros(k,k)

and

B=Q_{s}^{\top}A_{i}P_{s}

6 Set

i=1

P(2:k,1)=B(2:k,1)/R_{s}(1,1)

8 while $i<k-1$ do

i=i+1

P(i+1:k,i)=(B(i+1:k,i)-\sum\limits_{t=1}^{i-1}R_{s}(t,i)P(i+1:k,t)/R_{s}(i,i)

P(i,i+1:k)=-P(i+1:k,i)^{\top}

14 Set

Q_{i}=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}R_{s}^{\dagger}+Q_{s}P

15 Set

R_{i}=B-PR_{s}

Output:

Q=Q_{s}+Q_{i}\epsilon\in\mathbb{DR}^{m\times k}

is a dual matrix with orthogonal columns,

R=R_{s}+R_{i}\epsilon\in\mathbb{DR}^{k\times n}

is an upper trapezoidal dual matrix, and

P_{s}\in\mathbb{R}^{n\times n}

is a permutation matrix.

Algorithm 6 Randomized QR Decomposition with Column Pivoting in the standard part of Dual Matrices (RDQRCP)

4 Perturbation Analysis of the Q-factor

There has been extensive research [24, 26] on perturbation analysis for the thin QR decomposition of a matrix $A\in\mathbb{R}^{m\times n}$ with full column rank, where $A=QR,Q\in\mathbb{R}^{m\times n}$ and $R\in\mathbb{R}^{n\times n}$ . Suppose that $\Delta A\in\mathbb{R}^{m\times n}$ is a small perturbation matrix such that $A+\Delta A$ has full column rank, then $A+\Delta A$ has a unique QR factorization

A+\Delta A=(Q+\Delta Q)(R+\Delta R),

(4.1)

where $Q+\Delta Q$ has orthonormal columns and $R+\Delta R$ is upper triangular with positive diagonal entries. Perturbation analysis for the Q-factor in the QR factorization has received extensive attention in the literature [24, 26]. The infinitesimal part within dual numbers possesses highly distinctive properties

\epsilon^{2}=0,\qquad(\epsilon\neq 0).

(4.2)

Hence, if we view the decomposition of dual number matrices as the decomposition of matrices with perturbed terms, then what we obtain is an accurate first-order approximation. We substitute $A=A_{s}+\epsilon A_{i}\in\mathbb{DR}^{m\times n}$ for $A+\Delta A$ in (4.1), $Q=Q_{s}+\epsilon Q_{i}\in\mathbb{DR}^{m\times n}$ for $Q+\Delta Q$ , and $R=R_{s}+\epsilon R_{i}\in\mathbb{DR}^{n\times n}$ for $R+\Delta R$ . Hence, using dual numbers, we can obtain the explicitly accurate first-order perturbation [14] of the QR factorization,

A_{s}+\epsilon A_{i}=(Q_{s}+\epsilon Q_{i})(R_{s}+\epsilon R_{i}).

(4.3)

This allows us to view the QR decomposition of dual matrices from the perspective of perturbation analysis. For the three types of perturbation analyses previously mentioned and given conditions, it is evident that we can satisfy all of them. Stewart [24] initially provided perturbation bounds for the Q-factor in QR decomposition under componentwise computation,

\|Q_{i}\|_{F}\lesssim(1+\sqrt{2})\|A_{s}^{\dagger}\|_{2}\|A_{i}\|_{F}.

(4.4)

Subsequently, Sun [26] improved rigorous perturbation bounds for the Q- factor $Q_{i}$ when $A_{i}$ is given,

\|Q_{i}\|_{F}\lesssim\sqrt{2}\|A_{s}^{\dagger}\|_{2}\|A_{i}\|_{F}.

(4.5)

Now, we utilize the examples in [26] to test whether our tDQR algorithm satisfies the perturbation upper bounds as indicated in (4.5). Let

A_{s}=\begin{bmatrix}1&-2&1&2&3\\ 0&2&4&1&-5\\ 0&0&3&-1&2\\ 0&0&0&4&1\\ 0&0&0&0&5\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\end{bmatrix},\quad A_{i}=\tau\begin{bmatrix}0.2&-0.5&0.3&0.1&0.4\\ -0.1&0.4&0.1&-0.3&0.2\\ 0.5&0.7&-0.2&0.1&0.6\\ 0.3&-0.6&0.1&-0.1&0.2\\ 0.2&0.1&0.7&0.3&-0.4\\ 0.4&0.8&-0.2&0.1&0.3\\ 0.6&-0.1&-0.5&0.1&-0.2\\ 0.1&-0.3&0.2&0.6&0.7\end{bmatrix},

(4.6)

where $\tau>0$ . Table 1 illustrates the perturbation errors obtained through the MATLAB QR algorithm and our tDQR decomposition in Algorithm 3 for different small values of $\tau$ . Clearly, through numerical experiments, we can observe that the perturbation errors obtained using dual matrices are smaller and adhere to the perturbation bounds given by (4.5). This further supports the notion that dual matrix QR decomposition is an explicitly accurate first-order approximation with improved orthogonality.

Table 1: Perturbation Errors in tDQR Decomposition and QR Decomposition

$\tau$	1e-1	1e-2	1e-5	1e-8
$\\|A_{i}\\|_{F}$	0.24310492	0.02431049	2.43104916e-05	2.43104915e-08
$\\|\Delta Q\\|_{F}$	0.33558049	0.03162668	3.13819444e-05	3.13816983e-08
$\\|Q_{i}\\|_{F}$	0.31381698	0.03138170	3.13816980e-05	3.13816980e-08
$\sqrt{2}\\|A_{s}^{\dagger}\\|_{2}\\|A_{i}\\|_{F}$	1.04034046	0.10403405	1.04034046e-04	1.04034045e-07

5 Numerical Experiments

In this section, we compare the various proposed QR algorithms for the dual matrix and assess their respective suitability for different scenarios. Next, we utilize the QR decomposition of dual matrices to construct and prove the existence of dual Moore-Penrose generalized inverse (DMPGI). Numerical examples demonstrate that our approach achieves higher precision compared to previous methods. Subsequently, we verify the effectiveness of QR algorithms for the dual matrix in traveling wave detection for composite waves. Finally, we apply the proposed QR algorithm of the dual matrix to traveling wave detection in large-scale brain functional magnetic resonance imaging (fMRI) data in the field of neurosciences, significantly enhancing the efficiency relative to previous algorithms. The code was implemented using Matlab R2023b and all experiments were conducted on a personal computer equipped with an M2 pro Apple processor and 16GB of RAM, running macOS Ventura 13.3.

5.1 Numerical Comparison

In this section, we compare the runtime of the DQR, the DQRCP, the tDQ, and the tDQRCP decompositions of a dual matrix $A\in\mathbb{DR}^{m\times n}$ under various scenarios where $m>n,m=n$ , and $m<n$ , to assess their performance at different scales. We construct a GIID dual matrix $A$ as the synthetic data for evaluating and comparing. The assessment criteria comprise runtime measurements due to the absence of uniform norms, rendering computational error comparisons challenging. To ensure the reliability of our findings, each algorithm is executed 10 times, resulting in robust results.

Table 2: Comparison between DQR decomposition and DQRCP decomposition

	DQR	DQRCP
$500\times 500$	0.1267(1.1504e-4)	0.1344(3.5918e-5)
$500\times 1000$	0.1456(2.3226e-4)	0.1748(1.9633e-4)
$1000\times 500$	0.2727(0.0011)	0.2815(6.240e-4)
$2500\times 2500$	5.7008(0.1042)	5.7608(0.0054)
$2500\times 5000$	6.9206(0.0210)	8.0753(0.0622)
$5000\times 2500$	16.6889(0.4629)	17.4049(0.1718)
$4000\times 4000$	19.1934(1.4439)	22.3248(2.2537)
$4000\times 8000$	21.4972(0.6180)	28.0318(0.6010)
$8000\times 4000$	70.5938(0.9510)	73.6737(0.9375)
$5000\times 5000$	33.3352(1.4165)	38.3178(1.2820)
$5000\times 10000$	38.3174(2.1581)	51.5523(1.2879)
$10000\times 5000$	123.6001(18.5784)	133.7725(15.5367)

Table 3: Comparison between tDQR decomposition and tDQRCP decomposition

	tDQR	tDQRCP
$500\times 500$	0.1319(2.9287e-4)	0.1422(6.0144e-5)
$500\times 1000$	1.9193(0.0134)	2.9907(0.0430)
$1000\times 500$	0.1849(7.980e-4)	0.1808(3.636e-4)
$2500\times 2500$	5.7371(0.0094)	6.0572(0.0162)
$2500\times 5000$	18.1065(4.5056)	20.0397(1.7046)
$5000\times 2500$	8.3085(0.1454)	8.6551(0.1427)
$4000\times 4000$	20.4515(1.1058)	22.6906(1.5815)
$4000\times 8000$	62.7255(38.3143)	67.9553(20.5681)
$8000\times 4000$	28.2839(1.1321)	31.9453(0.9432)
$5000\times 5000$	34.1231(3.4861)	39.2982(1.8496)
$5000\times 10000$	109.0411(22.0643)	120.6198(20.0132)
$10000\times 5000$	48.1549(7.5332)	55.9513(6.7111)

Tables 2 and 3 present a time comparison for the DQR, the DQRCP, the tDQR, and the tDQRCP decompositions for dual matrices of varying dimensions, with variance values provided in parentheses. It is evident that QR decompositions with column pivoting consistently exhibit slower runtime than those without column pivoting. However, they offer enhanced numerical stability. When considering dual matrices of size $m\times n$ , we observe that for $m>n$ , both DQR decomposition and DQRCP decomposition outperform the tDQR decomposition and the tDQRCP decomposition in terms of speed. Conversely, when $m<n$ , the situation is reversed, thus confirming our hypothesis. When $m=n$ , it is noticeable that (3.1) lacks the presence of the term $(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}$ compared to (3.17). However, in cases where $m$ and $n$ are equal, and the dual matrix has full column rank, then $(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}=O_{m\times n}$ . Nonetheless, thin QR decomposition still computes this term, leading to increased runtime.

Next, we aim to construct a low-rank dual matrix $A\in\mathbb{DR}^{m\times n}$ , where $m>n$ . Let $L$ be a GIID $m\times r$ dual matrix, and $R$ be a GIID $r\times n$ dual matrix, where $r=round(n/10)$ . For the synthetic data $A=LR$ that we have constructed, we compare the efficiency of the low-rank decomposition using the tDQRCP decomposition and the RDQRCP decomposition, with the target rank $k$ set to 10. To ensure the reliability of our findings, each algorithm is executed 50 times, resulting in robust results.

Table 4: Comparison between tDQRCP decomposition and RDQRCP decomposition

	tDQRCP	RDQRCP
$1000\times 200$	0.0546(2.1961e-4)	0.0100(3.2969e-5)
$2000\times 400$	0.2088(0.0017)	0.0386(3.0649e-4)
$4000\times 1000$	1.4337(0.0114)	0.1557(2.1831e-4)
$8000\times 2000$	7.3579(0.0973)	0.6936(0.0058)

Table 4 presents a time comparison for the tDQRCP decomposition and the RDQRCP decomposition for dual matrices of varying dimensions, with variance values provided in parentheses. Clearly, in the low-rank problem, the RDQRCP decomposition exhibits faster speed and superior stability, with the potential to handle larger-scale data within memory constraints.

5.2 Dual Moore-Penrose Generalized Inverse

In previous research on dual matrices, many scholars have investigated the theory of dual Moore-Penrose generalized inverse (DMPGI) [1, 28, 29, 31]. We utilize the QR decomposition of dual matrices to compute DMPGI as an example. Similar to how Definition 2.1 defines the Moore-Penrose pseudoinverse, we can define the DMPGI accordingly.

Definition 5.1.

Given a dual matrix $A\in\mathbb{DR}^{m\times n}$ , its dual Moore-Penrose generalized inverse (DMPGI), denoted as $A^{\dagger}$ which is a dual matrix that satisfies the following four properties,

AA^{\dagger}A=A,\quad A^{\dagger}AA^{\dagger}=A^{\dagger},\quad(AA^{\dagger})^{\top}=AA^{\dagger},\quad(A^{\dagger}A)^{\top}=A^{\dagger}A.

(5.1)

In the following theorem, we establish the relationship between the existence of DMPGI and QR decomposition of the dual matrix.

Theorem 5.1.

Given a dual matrix $A\in\mathbb{DR}^{m\times n}$ where $A_{s}\in\mathbb{R}^{m\times n}$ is full column rank. Then the DMPGI $A^{\dagger}$ of $A$ exists, and

A^{\dagger}=R_{s}^{-1}Q_{s}^{\top}+(R_{s}^{-1}Q_{i}^{\top}-R_{s}^{-1}R_{i}R_{s}^{-1}Q_{s}^{\top})\epsilon

(5.2)

where $A=(Q_{s}+Q_{i}\epsilon)(R_{s}+R_{i}\epsilon)$ is the thin QR decomposition of $A$ .

Proof.

According to Theorem 3.2, since $A_{s}$ is full column rank, the thin QR decomposition of $A=QR=(Q_{s}+Q_{i}\epsilon)(R_{s}+R_{i}\epsilon)=Q_{s}R_{s}+(Q_{s}R_{i}+Q_{i}R_{s})\epsilon$ exists. Then the DMPGI of $A$ exists, taking the form $A^{\dagger}=R^{\dagger}Q^{\top}$ , satisfying Definition 5.1. Through rigorous computation, the matrix $R^{\dagger}=R_{s}^{-1}-R_{s}^{-1}R_{i}R_{s}^{-1}\epsilon$ can be derived, from which we can obtain

$\displaystyle A^{\dagger}$	$\displaystyle=R^{\dagger}Q^{\top}$	(5.3)
	$\displaystyle=(R_{s}^{-1}-R_{s}^{-1}R_{i}R_{s}^{-1}\epsilon)(Q_{s}^{\top}+Q_{i}^{\top}\epsilon)$
	$\displaystyle=R_{s}^{-1}Q_{s}^{\top}+(R_{s}^{-1}Q_{i}^{\top}-R_{s}^{-1}R_{i}R_{s}^{-1}Q_{s}^{\top})\epsilon.$

Based on the results provided by Theorem 5.1, we can compute the DMPGI using the QR decomposition of the dual matrix. Utilizing the examples in reference [16], let’s consider the dual matrix

A_{1}=\begin{bmatrix}1&3\\ 9&22\\ 4&4\end{bmatrix}+\epsilon\begin{bmatrix}4&0\\ 2&4\\ 4&1\end{bmatrix},\quad A_{2}=\begin{bmatrix}1&3&4\\ 9&22&4\end{bmatrix}+\epsilon\begin{bmatrix}4&0&1\\ 2&4&4\end{bmatrix}.

(5.4)

By the computation based on Theorem 5.1, we determine that the DMPGIs of $A_{1}$ and $A_{2}$ are respectively

A_{1}^{\dagger}=\begin{bmatrix}-0.0508&-0.0691&0.4182\\ 0.0276&0.0727&-0.1704\end{bmatrix}+\epsilon\begin{bmatrix}0.8221&-0.0350&-0.4596\\ -0.3493&0.0117&0.1675\end{bmatrix}

(5.5)

and

A_{2}^{\dagger}=\begin{bmatrix}-0.0349&0.0210\\ -0.0379&0.0438\\ 0.2872&-0.0381\\ \end{bmatrix}+\epsilon\begin{bmatrix}-0.0143&0.0002\\ -0.0350&-0.0011\\ -0.0071&-0.0107\\ \end{bmatrix}.

(5.6)

Remark 5.1.

Note that the DMPGI of $A_{1}$ calculated here differs from the result in reference [16]. Through simple numerical verification, we ascertain that our method provides higher accuracy. For the DMPGI of $A_{2}$ , our computation matches the result in reference [16], indicating that with the assistance of QR decomposition, we achieve higher precision in DMPGI. Our approach offers an effective, concise, and precise method for computing DMPGI.

5.3 Simulations of Standing and Traveling Waves

In the past, numerous theoretical explorations have paved the way for research into the application aspects of dual numbers. Eduard Study [25] introduced the concept of dual angles when determining the relative positions of two skew lines in three-dimensional space, where the angle is defined as the standard part and the distance as the infinitesimal part. Similarly, when dealing with raw time series data, we regard the standard part as the observed data and the infinitesimal part as its derivative or first-order difference, which naturally identifies standing and traveling waves in the simulations. Consider a spatiotemporal propagating pattern [8]

\textbf{x}(t)=2e^{\gamma t}[\cos{(\omega t)}\textbf{c}-\sin{(\omega t)}\textbf{d}],

(5.7)

where $\textbf{x}(t)$ represents a real vector of particle positions at the time $t$ , $\gamma$ , $\omega$ are real scalars representing the exponential decay rate and angular frequency, and c, d are real vectors representing two modes of oscillation. Overall, (5.7) demonstrates the circularity and continuity of propagation between c and d. When $\textbf{c}=\textbf{d}$ , (5.7) represents a ‘standing wave’ mode oscillation of rank-1. When $\textbf{c}\neq\textbf{d}$ , (5.7) represents a ‘traveling wave’ mode oscillation of rank-2 which behaves as a wavy motion. In previous research [8], wave propagation percentages of a wave were measured using an index, without distinguishing the exact traveling wave from the original wave. Wei, Ding, and Wei [34] introduced an effective method for traveling wave identification and extraction, including the determination of its propagation path. Now, we build upon their ideas and explore whether the QR decompositions of dual matrices proposed in this paper can be effectively utilized for traveling wave identification and whether they offer computational advantages. Here, we employ proper orthogonal decomposition to investigate wave behavior. Therefore, our focus lies in the $Q$ component of the QR decomposition for the dual matrix.

Refer to caption — Figure 1: The simulations of pure standing waves exhibit the center positions in the standard part of $Q$ and small values in the infinitesimal part. Simulations of pure traveling waves reveal paired similarities between the standard and infinitesimal parts of $Q$ . Leveraging the QR decomposition of the dual matrix and the specific properties of traveling waves, we accurately identify traveling waves within a combination wave composed of four standing waves and two traveling waves.

Now, we consider the spatiotemporal propagation patterns in (5.7) and construct a two-dimensional grid constructed as a $50\times 100$ array of pixels by Gaussian waves. We arrange $m$ particles in rows and $n$ time points in columns, allowing us to construct an $m\times n$ ensemble matrix $X_{s}$ . Here, we set c and d as two-dimensional Gaussian curves. By taking the derivative of $X_{s}$ with respect to time, we obtain $X_{i}$ . By treating $X_{s}$ as the standard part and $X_{i}$ as the infinitesimal part, we can construct the dual matrix $X$ . When $c=d$ , Figure 1(a) illustrates the standard part $Q_{s}$ and infinitesimal part $Q_{i}$ of the dual matrix $X$ after the QR decomposition. This is clearly a standing wave, with the center point at $(25,50)$ . When $c\neq d$ , Figure 1(b) displays the standard part $Q_{s}$ and infinitesimal part $Q_{i}$ of the dual matrix $X$ after QR decomposition, with center points at $(25,20)$ and $(25,80)$ . It’s noteworthy that signals $Q_{s}(:,1),Q_{i}(:,2),Q_{s}(:,2)$ , and $Q_{i}(:,1)$ exhibit paired similarities, where one is homologous and the other is heterologous. This property can be succinctly explained under the rank-2 truncated QR decomposition. First, due to the construction of $X$ , we have the range spaces $\textbf{Ran}(X_{i})\subseteq\textbf{Ran}(X_{s})$ , implying the existence of a matrix $M\in\mathbb{R}^{n\times n}$ such that $X_{s}=X_{i}M$ . Consequently, we have

(I_{m}-Q_{s}Q_{s}^{\top})A_{i}=(I_{m}-Q_{s}Q_{s}^{\top})A_{s}M=(I_{m}-Q_{s}Q_{s}^{\top})Q_{s}R_{s}M=O_{m\times n}.

(5.8)

For a $2\times 2$ skew-symmetric matrix $P=\begin{bmatrix}0&p_{12}\\ -p_{12}&0\end{bmatrix}$ , it is evident that

Q_{i}=Q_{s}P=\begin{bmatrix}Q_{s}(:,1)&Q_{s}(:,2)\end{bmatrix}\begin{bmatrix}0&p_{12}\\ -p_{12}&0\end{bmatrix}=\begin{bmatrix}-p_{12}Q_{s}(:,2)&p_{12}Q_{s}(:,1)\end{bmatrix}.

(5.9)

By utilizing (5.8) and (5.9) to analyze the $Q$ components of the QR decompositions proposed in this paper, it becomes apparent that paired similarities are an inevitable occurrence.

To validate the effectiveness of traveling wave identification, we construct a combination wave and attempt to separate different signals. We combine four standing waves and two traveling waves with different frequencies and weights, each with the same standard deviation $\sigma=1$ . Additionally, we introduce Gaussian noise with a peak value of $10^{-3}$ . Using the same methodology as before, we treat the data matrix of the combination wave as the standard part and its first-order derivative as the infinitesimal part to generate the dual matrix $Y$ . We perform a QR decomposition on $Y$ and utilize $Q$ to separate the six different waves. Figure 1(c) illustrates the standard and infinitesimal parts of the first four columns of $Q$ , corresponding to four weighted standing waves. The centers of these four standing waves (50, 50), (100, 100), (150, 70), and (170, 180) represent the wave peaks. Similarly, in the standard and infinitesimal parts of columns 5 to 8 of $Q$ , we can observe two pairs of traveling waves. In one pair, (50, 100) and (100, 50) represent the wave peaks, while in the other pair, (70, 150) and (120, 150) serve as the wave peaks.

In summary, our QR decomposition of the dual matrix provides an effective approach for traveling wave identification through proper orthogonal decomposition. Moreover, it is computationally straightforward.

5.4 Traveling Wave Identification in the Brain

Building upon prior research regarding the application of dual matrices in traveling wave identification for time series data in the brain, we hypothesize that performing a proper orthogonal decomposition on them yields superior results. Subsequently, we employ the previously proposed QR decomposition algorithm for dual matrices to test our hypothesis.

In our subsequent experiments, we made use of the Human Connectome Project (HCP) database [30], which comprises high-resolution 3D magnetic resonance scans of healthy adults aged between 22 to 35 years. This database includes task-state functional magnetic resonance imaging (tfMRI) images, a pivotal imaging modality in medical analysis. These tfMRI images were acquired while the subjects engaged in various tasks designed to stimulate cortical and subcortical networks. Each tfMRI scan encompasses two phases: one with right-to-left phase encoding and the other with left-to-right phase encoding. Achieving an in-plane field of view (FOV) rotation was accomplished by inverting both the readout (RO) and phase encoding (PE) gradient polarity ¹¹1https://www.humanconnectome.org/storage/app/media/documentation/s1200/HCP_S1200_Release_Reference_Manual.pdf.

Building upon the theoretical groundwork laid out in Section 5.3 regarding traveling wave identification by proper orthogonal decomposition, our present objective is to delve into the associations between this methodology and empirical brain responses within functional brain regions. Specifically, we concentrate on assessing the applicability of traveling wave identification within brain regions associated with typical language processing tasks encompassing all seven types, including semantic and phonological processing tasks, as initially developed by Binder et al. [4]. In this experimental setup, Binder et al. selected four blocks of a story comprehension task and four blocks of an arithmetic task. The story comprehension task involved the presentation of short auditory narratives consisting of 5-9 sentences, followed by binary-choice questions related to the central theme of the story. Conversely, the arithmetic task required participants to engage in addition and subtraction calculations during auditory experiments.

For our investigation, we considered a randomly chosen participant from the aforementioned experiment. The corresponding task-state functional magnetic resonance imaging (tfMRI) sequence data associated with the language processing task is stored as a 4-dimensional voxel-based image, encompassing $91\times 109\times 91$ three-dimensional spatial locations and 316 frames. The voxel size is set at 2mm, and the repetition time (TR) is established at 0.72 seconds. Data preprocessing assumes a pivotal role in maintaining substantial information while effectively compressing the image matrix. We initiated temporal filtering of the tfMRI signals, employing a Butterworth band-pass zero-phase filter in the frequency range of 0.01-0.1 Hz. This step aligns the signals with the conventional low-frequency range. Subsequently, we applied spatial smoothing, employing a 3-dimensional Gaussian low-pass filter with a half-maximum width of 5mm, involving convolution and exclusion of NANs. To ensure uniformity in voxel size across different frames, we implemented a Z-score transformation on the BOLD time series of all voxels. This transformation ensures a mean of zero and unit variance, maintaining consistency in the data.

Following this, our primary aim is to derive a dual matrix from the original brain sequence data matrix and affirm the coherence between areas displaying significant extracted traveling waves and the respective functional cortical activities in the brain. To achieve this, the initial three spatial dimensions are transformed into column vectors for each frame, and the fourth temporal dimension is sequentially appended. This operation culminates in the creation of a real matrix. This real matrix constitutes what we refer to as the standard component of the dual matrix $D$ . Similarly, the infinitesimal component is acquired by performing first-order temporal differentiation.

By utilizing the previously proposed RDQRCP in Algorithm 6 of the dual matrix $D$ , we easily obtain the desired proper orthogonal decomposition. Furthermore, compared to other researchers’ compact dual singular value decomposition (CDSVD) [34], our orthogonality is more pronounced. While CDSVD requires 14.661 seconds, we only need 5.660 seconds, providing a significant advantage in terms of processing time.

In Figure 2, the utilization of CDSVD [34] and the RDQRCP is presented in proper orthogonal decomposition, showcasing the results of inner products between the standard part and infinitesimal part. More specifically, within the QR decomposition for dual matrices, the inner product computed between $Q_{s}(:,16)$ and $Q_{i}(:,20)$ yields a substantial value of 0.840. Similarly, for $Q_{i}(:,20)$ and $Q_{s}(:,16)$ , the inner product registers -0.723. These calculations unveil the presence of a conspicuous traveling wave, characterized by the dominant principal components residing in the standard part. Subsequent to this discovery, the rank-2 matrix, originating from the second and third principal components, is meticulously rescaled to its original dimensions. This rescaled matrix is then subjected to spatial projection onto the brain’s anatomical structure using the ‘BrainNet’ toolkit within the MATLAB environment [35]. Through a sequential concatenation of these brain images, a dynamic video representation is created. This video provides a visually striking means of observing the temporal dynamics of brain activity across distinct regions, as depicted through the dynamic variations in color representation.

In Figure 3 (a), we can discern a substantial signal propagation. This module represents a language task in progress, where the wave originates from Brodmann area 39 (angular gyrus) and terminates at Brodmann area 40 (supramarginal gyrus) [22]. Both of these regions are intricately linked with Wernicke’s area, a pivotal node in language processing, particularly in the comprehension of speech. This specific cerebral region plays a critical role in the interpretation and generation of language, especially in speech comprehension. It enables us to adeptly grasp spoken language from others and articulate our own thoughts and emotions. Additionally, the prominence of this wave in the left hemisphere, as opposed to the right hemisphere, underscores the left hemisphere’s primary role in language processing.

Simultaneously, in Figure 3 (b), we also discern the propagation of the wave from Brodmann area 44 (pars opercularis) to Brodmann area 45 (Pars triangularis), both of which are constituents of Broca’s area [12]. This region plays a pivotal role in language processing, encompassing tasks related to the formulation and organization of language expression, as well as the processing of grammatical structures and syntax. Lesions or impairments in this area can lead to Broca’s aphasia, characterized by challenges in sentence production and reduced linguistic fluency, while language comprehension abilities remain largely intact. Additionally, it is implicated in language memory functions, particularly short-term memory processes.

To summarize, our findings have delineated two distinct wave patterns that correspond to well-established brain functional circuits, thus corroborating prior research in the field of cognitive neuroscience. In essence, this serves as a robust validation mechanism that reinforces the established conclusions within this domain.

Appendix A Proof of Theorem 3.2 and Theorem 3.3

The proof of Theorem 3.2.

Proof.

Suppose the thin QR decomposition of $A$ exists, then we have that

\begin{cases}A_{s}=Q_{s}R_{s}\\ A_{i}=Q_{s}R_{i}+Q_{i}R_{s}\end{cases}

(A.1)

through simple computation. By treating $R_{i}$ as $X$ and $Q_{i}$ as $Y$ , the given equation $A_{i}=Q_{s}R_{i}+Q_{i}R_{s}$ in (A.1) transforms into a Sylvester equation involving $Q_{i}$ and $R_{i}$ . Utilizing Theorem 2.1, if and only if

(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}(I_{n}-R_{s}^{\dagger}R_{s})=O_{m\times n},

(A.2)

which holds true identically by the invertibility nature of $R_{s}$ and the fact $R_{s}^{\dagger}=R_{s}^{-1}$ , one can derive the solution

\begin{cases}Q_{i}&=-(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}(-R_{s})^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\dagger})Z(-R_{s})(-R_{s})^{\dagger}\\ R_{i}&=Q_{s}^{\dagger}A_{i}+Q_{s}^{\dagger}Z(-R_{s})+(I_{n}-Q_{s}^{\dagger}Q_{s})W\end{cases}

(A.3)

where $Z\in\mathbb{R}^{m\times n}$ and $W\in\mathbb{R}^{n\times n}$ are arbitrary.

Considering the fact $Q_{s}^{\dagger}=Q_{s}^{\top}$ and $R_{s}^{\dagger}=R_{s}^{-1}$ , we have

$\displaystyle Q_{i}$	$\displaystyle=-(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}(-R_{s})^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\dagger})Z(-R_{s})(-R_{s})^{\dagger}$	(A.4)
	$\displaystyle=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}+Z-(I_{m}-Q_{s}Q_{s}^{\top})ZR_{s}R_{s}^{-1}$
	$\displaystyle=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}+Q_{s}Q_{s}^{\top}Z,$

and

	$\displaystyle R_{i}$	$\displaystyle=Q_{s}^{\dagger}A_{i}+Q_{s}^{\dagger}Z(-R_{s})+(I_{n}-Q_{s}^{\dagger}Q_{s})W$		(A.5)
		$\displaystyle=Q_{s}^{\top}A_{i}-Q_{s}^{\top}ZR_{s},$		(A.5)

where $Z\in\mathbb{R}^{m\times n}$ is arbitrary.

Since $Q_{s}\in\mathbb{R}^{m\times n}$ is a matrix with orthogonal columns, there exists a matrix $\widehat{Q}_{s}\in\mathbb{R}^{m\times(m-n)}$ with orthogonal columns, such that $[Q_{s},\widehat{Q}_{s}]$ is an orthogonal matrix. Thus, $Z=Q_{s}P+\widehat{Q}_{s}\widehat{P}$ , where $P\in\mathbb{R}^{n\times n}$ and $\widehat{P}\in\mathbb{R}^{(m-n)\times n}$ are arbitrary matrices. At this moment, we have $Q_{s}^{\top}Z=Q_{s}^{\top}(Q_{s}P+\widehat{Q}_{s}\widehat{P})=P$ . Then the expressions for $Q_{i}$ and $R_{i}$ are as follows,

\begin{cases}Q_{i}&=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}+Q_{s}P\\ R_{i}&=Q_{s}^{\top}A_{i}-PR_{s},\end{cases}

(A.6)

where $P\in\mathbb{R}^{n\times n}$ is an arbitrary matrix.

Just as we discussed in Theorem 3.1, we aim to exploit the orthogonal column properties of $Q$ and the upper triangular nature of $R$ to ascertain the values of $P$ . In this regard, we first contemplate the scenario where $Q$ is a dual matrix with orthogonal columns. Thus, we have

Q_{s}^{\top}Q_{i}+Q_{i}^{\top}Q_{s}=O_{n}.

(A.7)

By substituting the expression for $Q_{i}$ from (A.6), we can obtain

$\displaystyle O_{n}$	$\displaystyle=Q_{s}^{\top}[(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}+Q_{s}P]+[(I_{m}-Q_{s}Q_{s}^{\top})A_{i}R_{s}^{-1}+Q_{s}P]^{\top}Q_{s}$	(A.8)
	$\displaystyle=Q_{s}^{\top}Q_{s}P+(Q_{s}P)^{\top}Q_{s}$
	$\displaystyle=P+P^{\top}$

with the fact $Q_{s}^{\top}(I_{m}-Q_{s}Q_{s}^{\top})=O_{n}$ .

We deduce that $P\in\mathbb{R}^{n\times n}$ is a skew-symmetric matrix from (A.8). Next, we leverage the upper triangular property of $R$ to ascertain the elements within matrix $P$ . Considering the expression of $R_{i}$ in (A.6), we have

PR_{s}=Q_{s}^{\top}A_{i}-R_{i}.

(A.9)

We denote $B=Q_{s}^{\top}A_{i}\in\mathbb{R}^{n\times n}$ and $B=[b_{1},b_{2},\ldots,b_{n}],R_{i}=[r_{i1},r_{i2},\ldots,r_{in}],R_{s}=[r_{s1},r_{s2},\ldots,r_{sn}]$ where $b_{k},r_{ik},r_{sk}\in\mathbb{R}^{n},k=1,2,\ldots,n$ . By partitioning the matrices column-wise and given that $R_{s}$ is of full rank, (A.9) is transformed into the following system of $n$ equations,

\begin{cases}Pr_{s1}&=b_{1}-r_{i1},\\ Pr_{s2}&=b_{2}-r_{i2},\\ &\ldots\\ Pr_{sn}&=b_{n}-r_{in}.\\ \end{cases}

(A.10)

Exploiting that matrix $R_{s}$ is upper triangular, along with $R_{s}=(r_{s_{ij}})\in\mathbb{R}^{n\times n}$ and $P=[p_{1},p_{2},\ldots,p_{n}]$ where $p_{k}\in\mathbb{R}^{n},k=1,2,\ldots n$ , the system (A.10) can be transmuted into

\begin{cases}r_{s_{11}}p_{1}&=b_{1}-r_{i1},\\ r_{s_{12}}p_{1}+r_{s_{22}}p_{2}&=b_{2}-r_{i2},\\ &\ldots\\ r_{s_{1n}}p_{1}+r_{s_{2n}}p_{2}+\ldots+r_{s_{nn}}p_{n}&=b_{n}-r_{in}.\end{cases}

(A.11)

Recognizing the upper triangular structure of $R_{i}$ , it follows that for $r_{ik},k=1,2,\ldots,(n-1)$ , the lower $n-k$ entries $r_{ik}(k+1:n)$ are all zeros. Consequently, $p_{k}$ also requires computation solely for the lower $n-k$ entries $p_{k}(k+1:n)$ , thereby ensuring the maintenance of $P$ ’s skew-symmetric property. Hence, the system of $(n-1)$ equations has the solutions

\begin{cases}p_{1}(2:n)&=b_{1}(2:n)/r_{s_{11}},\\ p_{2}(3:n)&=(b_{2}(3:n)-r_{s_{12}}p_{1}(3:n))/r_{s_{22}},\\ &\ldots\\ p_{n-1}(n)&=(b_{n-1}(n)-\sum\limits_{t=1}^{n-2}r_{s_{t(n-1)}}p_{t}(n))/r_{s_{(n-1)(n-1)}}.\end{cases}

(A.12)

Thus, we complete the proof.

The proof of Theorem 3.3.

Proof.

Suppose randomized QR decomposition with column pivoting in the standard part of $A$ exists, then we have that

\begin{cases}A_{s}P_{s}=Q_{s}R_{s}\\ A_{i}P_{s}=Q_{s}R_{i}+Q_{i}R_{s}\end{cases}

(A.13)

through computation directly. By treating $R_{i}$ as $X$ and $Q_{i}$ as $Y$ , the given equation $A_{i}P_{s}=Q_{s}R_{i}+Q_{i}R_{s}$ in (A.13) transforms into a Sylvester equation involving $Q_{i}$ and $R_{i}$ . Utilizing Theorem 2.1 and the fact $Q_{s}^{\dagger}=Q_{s}^{\top}$ , if and only if

(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}(I_{n}-R_{s}^{\dagger}R_{s})=O_{m\times n}

(A.14)

holds true, one can derive the solution

\begin{cases}Q_{i}&=-(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}P_{s}(-R_{s})^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\dagger})Z(-R_{s})(-R_{s})^{\dagger}\\ R_{i}&=Q_{s}^{\dagger}A_{i}P_{s}+Q_{s}^{\dagger}Z(-R_{s})+(I_{k}-Q_{s}^{\dagger}Q_{s})W\end{cases}

(A.15)

where $Z\in\mathbb{R}^{m\times k}$ and $W\in\mathbb{R}^{k\times n}$ are arbitrary.

Considering the fact $Q_{s}^{\dagger}=Q_{s}^{\top}$ and $R_{s}R_{s}^{\dagger}=I_{k}$ , we have

$\displaystyle Q_{i}$	$\displaystyle=-(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}P_{s}(-R_{s})^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\dagger})Z(-R_{s})(-R_{s})^{\dagger}$	(A.16)
	$\displaystyle=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}R_{s}^{\dagger}+Z-(I_{m}-Q_{s}Q_{s}^{\top})ZR_{s}R_{s}^{\dagger}$
	$\displaystyle=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}R_{s}^{\dagger}+Q_{s}Q_{s}^{\top}Z,$

and

	$\displaystyle R_{i}$	$\displaystyle=Q_{s}^{\dagger}A_{i}P_{s}+Q_{s}^{\dagger}Z(-R_{s})+(I_{n}-Q_{s}^{\dagger}Q_{s})W$		(A.17)
		$\displaystyle=Q_{s}^{\top}A_{i}P_{s}-Q_{s}^{\top}ZR_{s},$		(A.17)

where $Z\in\mathbb{R}^{m\times k}$ is arbitrary.

Since $Q_{s}\in\mathbb{R}^{m\times k}$ is a matrix with orthogonal columns, there exists a matrix $\widehat{Q}_{s}\in\mathbb{R}^{m\times(m-k)}$ with orthogonal columns, such that $[Q_{s},\widehat{Q}_{s}]$ is an orthogonal matrix. Thus, $Z=Q_{s}P+\widehat{Q}_{s}\widehat{P}$ , where $P\in\mathbb{R}^{k\times k}$ and $\widehat{P}\in\mathbb{R}^{(m-k)\times k}$ are arbitrary matrices. At this moment, we have $Q_{s}^{\top}Z=Q_{s}^{\top}(Q_{s}P+\widehat{Q}_{s}\widehat{P})=P$ .Then the expressions for $Q_{i}$ and $R_{i}$ are as follows,

\begin{cases}Q_{i}&=(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}R_{s}^{\dagger}+Q_{s}P\\ R_{i}&=Q_{s}^{\top}A_{i}P_{s}-PR_{s},\end{cases}

(A.18)

where $P\in\mathbb{R}^{k\times k}$ is an arbitrary matrix.

First, we contemplate the scenario where $Q$ is a dual matrix with orthogonal columns. Thus, we have

Q_{s}^{\top}Q_{i}+Q_{i}^{\top}Q_{s}=O_{k}.

(A.19)

By substituting the expression for $Q_{i}$ from (A.18), we can obtain

$\displaystyle O_{k}$	$\displaystyle=Q_{s}^{\top}[(I_{m}-Q_{s}Q_{s}^{\top})A_{i}P_{s}R_{s}^{\dagger}+Q_{s}P]+[(I_{m}-Q_{s}Q_{s}^{\dagger})A_{i}P_{s}R_{s}^{\dagger}+Q_{s}P]^{\top}Q_{s}$	(A.20)
	$\displaystyle=Q_{s}^{\top}Q_{s}P+(Q_{s}P)^{\top}Q_{s}$
	$\displaystyle=P+P^{\top}$

with the fact $Q_{s}^{\top}(I_{m}-Q_{s}Q_{s}^{\top})=O_{k}$ .

We deduce that $P\in\mathbb{R}^{k\times k}$ is a skew-symmetric matrix from (A.20). Next, we leverage the upper trapezoidal property of $R$ to ascertain the elements within matrix $P$ . Considering the expression of $R_{i}$ in (A.18), we have

PR_{s}=Q_{s}^{\top}A_{i}P_{s}-R_{i}.

(A.21)

We denote $B=Q_{s}^{\top}A_{i}P_{s}\in\mathbb{R}^{k\times n}$ and $B=[b_{1},b_{2},\ldots,b_{n}],R_{i}=[r_{i1},r_{i2},\ldots,r_{in}],R_{s}=[r_{s1},r_{s2},\ldots,r_{sn}]$ where $b_{t},r_{it},r_{st}\in\mathbb{R}^{k},t=1,2,\ldots,n$ . By partitioning the matrices column-wise and given that $R_{s}$ is of full rank, (A.21) is transformed into the following system of $n$ equations,

\begin{cases}Pr_{s1}&=b_{1}-r_{i1},\\ Pr_{s2}&=b_{2}-r_{i2},\\ &\ldots\\ Pr_{sn}&=b_{n}-r_{in}.\\ \end{cases}

(A.22)

Exploiting that matrix $R_{s}$ is upper trapezoidal, along with $R_{s}=(r_{s_{ij}})\in\mathbb{R}^{k\times n}$ and $P=[p_{1},p_{2},\ldots,p_{k}]$ where $p_{t}\in\mathbb{R}^{k},t=1,2,\ldots n$ , the overdetermined system (A.22) can be transmuted into

\begin{cases}r_{s_{11}}p_{1}&=b_{1}-r_{i1},\\ r_{s_{12}}p_{1}+r_{s_{22}}p_{2}&=b_{2}-r_{i2},\\ &\ldots\\ r_{s_{1n}}p_{1}+r_{s_{2n}}p_{2}+\ldots+r_{s_{kn}}p_{k}&=b_{n}-r_{in}.\end{cases}

(A.23)

Recognizing the upper trapezoidal structure of $R_{i}$ , it follows that for $r_{it},t=1,2,\ldots,(k-1)$ , the lower $k-t$ entries $r_{it}(t+1:k)$ are all zeros. Consequently, $p_{t}$ also requires computation solely for the lower $k-t$ entries $p_{t}(t+1:k)$ , thereby ensuring the maintenance of $P$ ’s skew-symmetric property. Moreover, the consistency of an over-determined system (A.22) can be uniquely determined by $R_{i}(:,k+1:n)$ . Hence, the over-determined system of $(k-1)$ equations has the solutions

\begin{cases}p_{1}(2:k)&=b_{1}(2:k)/r_{s_{11}},\\ p_{2}(3:k)&=(b_{2}(3:k)-r_{s_{12}}p_{1}(3:k))/r_{s_{22}},\\ &\ldots\\ p_{k-1}(k)&=(b_{k-1}(k)-\sum\limits_{t=1}^{k-2}r_{s_{t(k-1)}}p_{t}(k))/r_{s_{(k-1)(k-1)}}.\end{cases}

(A.24)

Thus, we complete the proof.

Data Availability

The numerical data shown in the figures are available upon request, and so is the code to run the experiment and generate the figures.

Declarations

Conflicts of Interest. The authors declare no competing interests.

References

[1] J. Angeles. The dual generalized inverses and their applications in kinematic synthesis. In Latest advances in robot kinematics, pages 1–10. Springer, 2012.
[2] J. K. Baksalary and R. Kala. The Matrix Equation ${AX-YB=C}$ . Linear Algebra and its Applications, 25:41–43, 1979.
[3] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind. Automatic differentiation in machine learning: a survey. Journal of Marchine Learning Research, 18:1–43, 2018.
[4] J. R. Binder, W. L. Gross, J. B. Allendorfer, L. Bonilha, J. Chapin, J. C. Edwards, T. J. Grabowski, J. T. Langfitt, D. W. Loring, M. J. Lowe, et al. Mapping anterior temporal lobe language areas with fMRI: a multicenter normative study. Neuroimage, 54(2):1465–1475, 2011.
[5] M. Clifford. Preliminary sketch of biquaternions. Proceedings of the London Mathematical Society, 1(1):381–395, 1871.
[6] C. Cui, H. Wang, and Y. Wei. Perturbations of moore-penrose inverse and dual moore-penrose generalized inverse. Journal of Applied Mathematics and Computing, pages 1–24, 2023.
[7] J. A. Duersch and M. Gu. Randomized QR with column pivoting. SIAM Journal on Scientific Computing, 39(4):C263–C291, 2017.
[8] B. Feeny. A complex orthogonal decomposition for wave motion analysis. Journal of Sound and Vibration, 310(1-2):77–90, 2008.
[9] J. A. Fike and J. J. Alonso. Automatic differentiation through the use of hyper-dual numbers for second derivatives. In Recent Advances in Algorithmic Differentiation, pages 163–173. Springer, 2012.
[10] M. Fliess and C. Join. Model-free control. International Journal of Control, 86(12):2228–2252, 2013.
[11] B. Fornberg. Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of computation, 51(184):699–706, 1988.
[12] A. D. Friederici, N. Chomsky, R. C. Berwick, A. Moro, and J. J. Bolhuis. Language, mind and brain. Nature human behaviour, 1(10):713–722, 2017.
[13] G. H. Golub and C. F. Van Loan. Matrix Computations. JHU Press, 4th edition, 2013.
[14] A. Greenbaum, R.-C. Li, and M. L. Overton. First-order perturbation theory for eigenvalues and eigenvectors. SIAM Rev., 62(2):463–482, 2020.
[15] R. Gutin. Generalizations of singular value decomposition to dual-numbered matrices. Linear and Multilinear Algebra, 70(20):5107–5114, 2022.
[16] E. Pennestrì and R. Stefanelli. Linear algebra and numerical algorithms using dual numbers. Multibody System Dynamics, 18:323–344, 2007.
[17] E. Pennestrì, P. Valentini, and D. De Falco. The Moore–Penrose dual generalized inverse matrix with application to kinematic synthesis of spatial linkages. Journal of Mechanical Design, 140(10):102303, 2018.
[18] R. Peón, O. Carvente, C. A. Cruz-Villar, M. Zambrano-Arjona, and F. Peñuñuri. Dual numbers for algorithmic differentiation. Ingeniería, 23(3):71–81, 2019.
[19] L. Qi, D. M. Alexander, Z. Chen, C. Ling, and Z. Luo. Low rank approximation of dual complex matrices. arXiv:2201.12781, 2022.
[20] L. Qi and C. Cui. Eigenvalues and Jordan forms of dual complex matrices. Communications on Applied Mathematics and Computation, pages 1–17, 2023.
[21] L. Qi and Z. Luo. Eigenvalues and singular values of dual quaternion matrices. Pacific Journal of Optimization, 19(2):257–272, 2023.
[22] Y. Sakurai. Brodmann areas 39 and 40: human parietal association area and higher cortical function. Brain and nerve, 69(4):461–469, 2017.
[23] J. Sola. Quaternion kinematics for the error-state kalman filter. arXiv:1711.02508, 2017.
[24] G. Stewart. On the perturbation of LU, Cholesky, and QR factorizations. SIAM Journal on Matrix Analysis and Applications, 14(4):1141–1145, 1993.
[25] E. Study. Geometrie der Dynamen. Druck und verlag von BG Teubner, 1903.
[26] J.-G. Sun. On perturbation bounds for the QR factorization. Linear Algebra and its Applications, 215:95–111, 1995.
[27] F. E. Udwadia. Dual generalized inverses and their use in solving systems of linear dual equations. Mechanism and Machine Theory, 156:104158, 2021.
[28] F. E. Udwadia. When does a dual matrix have a dual generalized inverse? Symmetry, 13(8):1386, 2021.
[29] F. E. Udwadia, E. Pennestri, and D. de Falco. Do all dual matrices have dual Moore–Penrose generalized inverses? Mechanism and Machine Theory, 151:103878, 2020.
[30] D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. Behrens, E. Yacoub, K. Ugurbil, W.-M. H. Consortium, et al. The Wu-Minn human connectome project: an overview. Neuroimage, 80:62–79, 2013.
[31] H. Wang. Characterizations and properties of the MPDGI and DMPGI. Mechanism and Machine Theory, 158:104212, 2021.
[32] H. Wang, C. Cui, and X. Liu. Dual $r$ -rank decomposition and its applications. Computational and Applied Mathematics, 42(8):349, 2023.
[33] J. WB and J. LINDENSTRAUSS. Extensions of Lipschitz mappings into a Hilbert space. In Conference in Modern Analysis and Probability, pages 189–206, 1984.
[34] T. Wei, W. Ding, and Y. Wei. Singular value decomposition of dual matrices and its application to traveling wave identification in the brain. SIAM Journal on Matrix Analysis and Applications, 45(1):634–660, 2024.
[35] M. Xia, J. Wang, and Y. He. Brainnet viewer: a network visualization tool for human brain connectomics. PloS one, 8(7):e68910, 2013.
[36] Z. Xiong, Y. Wei, R. Xu, and Y. Xu. Low-rank traffic matrix completion with marginal information. Journal of Computational and Applied Mathematics, 410:114219, 2022.
[37] R. Xu, T. Wei, Y. Wei, and H. Yan. UTV decomposition of dual matrices and its applications. Computational and Applied Mathematics, 43(1):41, 2024.
[38] K. Yang and C. Shahabi. A pca-based similarity measure for multivariate time series. In Proceedings of the 2nd ACM international workshop on Multimedia databases, pages 65–74, 2004.