Hadamard matrices are n × n matrices H with elements from {-1, +1} such that H • H t = n • I n , where t denotes matrix transposition and I n is the n × n identity matrix. These matrices are considered trivial when their order is either n = 1 or n = 2. Beyond these two trivial orders, a well-known necessary condition for the existence of Hadamard matrices is that n ≡ 0 (mod 4). The sufficiency of this condition is the famous Hadamard conjecture, still open after more than 130 years [1,2]. Despite several dozens of known constructions for Hadamard matrices, the Hadamard conjecture remains open, mainly because each construction either covers a rather sparse set of orders or fails for several orders. There is a construction that furnishes a promising avenue towards a more structured form of the Hadamard conjecture, i.e., that conjecturally does not fail for any order which is a multiple of four [3]. This construction is based on a combinatorial object called Legendre pairs. Legendre pairs were introduced in [4], where the authors show how to apply the DFT matrix to search for these combinatorial objects. The importance of Legendre pairs is that they can be used to build Hadamard matrices and the conjecture that Legendre pairs of every odd length exist, implies the Hadamard conjecture, albeit via the very specialized structure of Legendre pairs [3]. From the computational point of view, the smallest odd order for which Legendre pairs are not known to exist is 115 [5]. In order to tackle the smallest open case, i.e., ℓ = 115 with our new algorithm, one would first have to further carefully integrate and incorporate the concept of Djokovic-Kotsireas compression [6] into the current implementation, as well as further optimize the current code. This is planned to be the objective of future work.
In this paper, we explore a new perspective to search for Legendre pairs, i.e., using matrix analysis and linear algebra. We start by exploring the structure of the Legendre pair matrix equation. Structural properties of matrices can be explored to solve many problems in applied sciences and engineering efficiently. These matrices have been studied individually in mathematics, engineering, and computer science disciplines for quite a long time. Within the last two decades, the theories of structured matrices got more attention as they are positioned to bridge diverse disciplines [7,8,9,10,11,12]. More specifically, recent work has been done to show that the theories of matrix analysis could be utilized to solve problems in narrowband multi-beam beamforming [13], wideband multi-beam beamforming [14], mutual coupling effects in antenna arrays [15], image processing [16,17], and reduction and minimization of circuit complexity and power consumption [18,13,19,20,14,21] through the derivation of FFT-like algorithms. Challenges remain to develop algorithms for solving problems in diverse disciplines while offering low-complexity, accuracy, and stability at the same time [9,22]. However, if we could solve linear systems using fast and recursive algorithms, we could obtain numerically stable algorithms with the same order of complexity [23,24,25,26]. Thus, we exploit the structure of the Legendre pair matrix equation to obtain a low-complexity and exact algorithm for searching for these objects.
The FFT is an algorithm that is used to compute the DFT matrix and its inverse efficiently.
The FFT algorithm can be used to reduce the brute-force calculation of the n × n DFT matrix by an n × 1 vector from O(n 2 ) to O(n log(n)), where n = 2 t (t ≥ 1) [27,28,29]. The structure and properties of the DFT matrix yield this complexity reduction. An n × n circulant is another structured matrix and can be defined using O(n) as opposed to O(n 2 ) parameters. More specifically, circulant matrices are completely determined using the first row (or column) by a cyclic shift of the entries to the right (or down). Furthermore, it is known that the circulant matrices can be diagonalized using the DFT and its conjugate transpose and hence computed by using the FFT algorithms [30,31,32,33]. In this paper, we introduce and incorporate the FFT technique to search for Legendre pairs of odd length ℓ, efficiently. Thus, we obtain visible speed-ups, especially for odd lengths > 50 while reducing the arithmetic complexity from O(ℓ 2 ) to O(n log(n)), where n > ℓ is the closest power of 2 to ℓ.
The organization along with the contribution of the paper is as follows. Section 2 presents preliminary results illustrating the use of structured matrices to study the Legendre pair matrix equation and its correlation with the established properties of the two circulant core (2cc) Hadamard matrix. Section 3 presents an algorithm with low arithmetic complexity for computing the product of the DFT matrix with each sequence of the Legendre pair. This algorithm allows for an efficient search for Legendre pairs while minimizing both arithmetic and time complexities. Section 4 presents numerical results demonstrating that the sparse factorization discussed in Section 3 yields an FFT-like algorithm. This, in turn, is utilized in a C implementation, resulting in the algorithm with a low time complexity for searching for Legendre pairs. This algorithm overreaches the time complexity of the existing combinatorial algorithms in searching for Legendre pairs. Finally, Section 5 concludes the paper.
This section begins by reviewing the definition of the Legendre pair matrix equation and discussing its basic properties. In the sequel, the Legendre pair matrix equation will be presented as a structured matrix. Legendre pairs are defined via the concept of the periodic autocorrelation function (PAF). Thus, we start the discussion with the definition of the PAF followed by the Legendre pair sequence as stated below.
where k + s is taken mod n.
Each PAF A (s) value is the inner product of the sequence L A and its right cyclic shift by s. For technical reasons and w.l.o.g. (see [34]) we impose the additional normalization conditions
where a i and b i , for i = 1, 2, . . ., ℓ denote the elements of L A , L B respectively.
When (L A , L B ) forms a Legendre pair of order ℓ and A, B are two ℓ × ℓ circulant matrices with their first rows as L A , L B respectively, a 2cc Hadamard matrix can be constructed as:
where + stands for +1 and -stands for -1. Let us recall that the determinant of the Hadamard matrix (2), is given by
In the sequel, we show a relationship between det(H 2ℓ+2 ) and the Legendre pair matrix equation, which will be defined next, to reassure the Property 3 based on the structured matrices. It is important to note that an even-order Legendre pair does not exist. Had there been one, we would have had a Hadamard matrix of order 2ℓ + 2, which is congruent to 2 modulo 4. Nevertheless, such a matrix does not exist.
This section starts with the Legendre pair matrix equation, followed by reestablishing the connection between the Legendre pair matrix equation and the 2cc Hadamard matrix via block Gauss-Jordan form and structured matrices.
Lemma 2.3. (Legendre pair matrix equation in [35]) Let (L A , L B ) be a Legendre pair of (odd) order ℓ. Let A, B be two ℓ × ℓ circulant matrices, whose first rows are L A , L B respectively.
Then the following matrix equation is satisfied:
where I ℓ denotes the ℓ × ℓ identity matrix and J ℓ denotes the ℓ × ℓ all-ones matrix. Moreover, the equation ( 4) is called the Legendre pair matrix equation.
We denote the right-hand-side matrix of the Legendre pair matrix equation by M (ℓ) = (2ℓ +
2)I ℓ -2J ℓ , and remark that M (ℓ) is a circulant matrix with first row equal to r = [2ℓ, -2, . . . -2
, where i, j = 1, 2, . . ., ℓ, -2 if i ̸ = j, where i, j = 1, 2, . . ., ℓ.
(
In the subsequent proposition, we show a relationship linking the determinants of 2cc
Hadamard matrices denoted as H 2ℓ+2 in equation (2), with the Legendre pair matrix equation (4), using block Gauss-Jordan form of the 2cc Hadamard matrix. Proposition 2.4. Let H 2ℓ+2 = H 11 H 12 H t 12 H 22 and A, B be two ℓ × ℓ circulant matrices satisfying Equation (4), where
, and
Proof. One can decompose the matrix H 22 into lower and upper block triangular matrices using the Schur complement of the block A t of the matrix H 22 s.t.
where A is the non-singular circulant matrix. Thus, we get
Similarly, we could decompose H 2ℓ+2 using the nonsingular block circulant matrix H 22 s.t.
and obtain det(
). On the other hand, we can continue the block lower and upper decomposition of H 22 in order to decompose the matrix H 22 using a block Gauss-Jordan form of H 22 s.t.
Since A is non-singular and H 22 admits the Gauss-Jordan elimination, we could obtain the inverse of H 22 s.t.
Thus, from the above, using the circulant matrices A and B generated by L A and L B , respectively, and followed by the block matrix multiplication, we get
Since A and B are circulant matrices derived from the Legendre pair L A and L B , respectively, with the normalization condition (1), the sum of each column in the circulant matrices (A t ) -1 B t and B(A t ) -1 equals 1. As a result, the entries in the first and last matrices in the factorization (8) will be simplified. Hence, from the above followed by the matrix multiplication, we get
From the above and using the matrix H 11 , we get
Thus, by ( 7) and ( 9), we get
Since A and B are non-singular circulant matrices, they are commutative, so we have
Thus, we get
Hence, from (4), (10), and (11), we get the result.
In the following proposition, we use the structure of the matrix M (ℓ) identified via the Legendre pair matrix equation ( 4) to obtain the determinant of the Legendre pair matrix equation, for any odd length ℓ.
Proposition 2.5. Let M (ℓ) be defined as in (5). Then,
Proof. The matrix M (ℓ) is circulant. Thus, by using the well-known decomposition of the circulant matrices via the DFT matrices [36,30,9,37,14], M (ℓ) can be expressed via
where
and F * ℓ is the conjugate transpose of F ℓ . Furthermore, the diagonal elements of D ℓ are eigenvalues of the circulant matrix M (ℓ) . Let these eigenvalues be represented as λ j 's associated with M (ℓ) , then these can be explicitly expressed as
where j = 0, 1, . . ., ℓ -1. Thus, the above equation is simplified into
Since F ℓ and F * ℓ are unitary matrices and following equations ( 13) and ( 15), we get the result.
We note that Propositions 2.4 and 2.5 confirm that the determinant of the Hadamard matrix can be calculated using equation (3).
In this section, we continue utilizing the structure of the Legendre pair matrix equation to obtain eigenvalues, eigenvectors, characteristic polynomial, spectral radius, and a bound for the spectrum of the matrices A and B. These properties will be utilized as auxiliary results for the facts on PSD values in Section 3.
Since the Legendre pair matrix equation can be seen as a circulant matrix M (ℓ) and the matrix M (ℓ) admits the decomposition (13), the normalized eigenvectors of the circulant matrix M (ℓ) are the Fourier modes given via v j = 1 √ ℓ (1, ω j , ω 2j , . . ., ω (ℓ-1)j ), for j = 0, 1, . . ., ℓ -1. Moreover, the eigenvalues of the Legendre pair matrix equation are given via equation (15). Thus, the matrix M (ℓ) has the eigenvalues 2 (with multiplicity 1) and 2ℓ + 2 (with multiplicity ℓ -1). This leads to the following immediate result based on the characteristic polynomial of the matrix M (ℓ) . Corollary 2.6. Let M (ℓ) be defined as in (5). Then its characteristic polynomial, denoted by
Proof. This immediately follows from equations ( 13) and (15).
Since the eigenvalues of the Legendre pair matrix equation are known, we can state the spectral radius as follows.
Corollary 2.7. Let M (ℓ) be defined as in (5). Then its spectral radius is given via
Proof. This is a direct consequence of the eigenvalues of M (ℓ) listed in (15).
Since the eigenvalues and spectral radius of the Legendre pair matrix equation are stated, we could obtain a bound for the spectrum of eigenvalues based on the Gershgorin disc. Before presenting the next Corollary, let us introduce the Gershgorin disc as follows [38,39].
Definition 2.8. (Gershgorin disc) Let A be a complex n × n matrix with entries a ij . For i ∈ {1, 2, . . ., n}, let R i be the sum of the absolute values of the non-diagonal entries in the i th row
In the following, we use the notation D(a, R) ⊆ C to denote a closed disc with center a ∈ C and radius R > 0.
Corollary 2.9. Let M (ℓ) be defined as in (5). Then its eigenvalues lie within the Gershgorin circle
Proof. This follows immediately from Propositions 2.4 and 2.5, and the Gershgorin circle theorem [38,39]. We utilize the following example to illustrate the above Corollary.
Example 2.10. Consider the Legendre pair of order ℓ = 45, given by the sequences L A and L B :
In Figure 1 Continuing from Corollary 2.9, we obtain the following result based on the eigenvalues of the
where
. Furthermore, the matrices D A and D B consist of eigenvalues of the matrices A and B, respectively. Let us denote the eigenvalues of A and B by (λ A ) j and (λ B ) j , respectively, for j = 0, 1, . . ., ℓ -1. Thus, using equation ( 18) allows us to express the entries of D A and D B by:
Thus, when j = 0, the above equations simplify to
Therefore, based on property (1), i.e., the normalization property yielding the sums in the RHS of the above equation to 1, we get that 1 is an eigenvalue of A and also B.
In this section, we first introduce the DFT matrix followed by the FFT and FFTW algorithms. Finally, we obtain bounds for the PSD values, which can be useful in further pruning the search space of Legendre pairs.
The FFT is the most frequently used algorithm in digital signal processing [40,41,42], and is used to efficiently compute the DFT and its inverse. The DFT of a sequence of n inputs x = {x k } n-1 k=0 is a sequence of n outputs y = {y k } n-1 k=0 defined via y = F n x, where
is the n × n normalized DFT matrix and ω n = e -2πi n , which is the superset of ω in (13) as ω defines only the primitive odd, i.e., ℓ th , root of unity. When n = 2 t (t ≥ 1), the DFT matrix can be computed using the FFT algorithm [27,28] because the DFT matrix admits the decomposition
where The difference in complexity of the FFT algorithms from 1965 [27] and the lowest best-known FFT [29] is about 25%, and yet there is a gap between mathematical theories and implementation of the highly optimized FFT packages [43,44]. Authors in these books implemented FFT in practice by creating the FFTW library, which is a widely used free software library that computes the DFT and its various special cases. Here, we will adopt the FFTW to compute the PSD of the candidate Legendre pairs with the reduction of arithmetic complexity from O(ℓ 2 ) to O(n log(n)), and also a drastic time complexity reduction in Section 4.
In this section, we will propose a low arithmetic complexity algorithm for the computation of the set of PSD values, as an alternative to the conventional brute-force approach involving the product of the DFT matrix with each of the two sequences L A and L B of the Legendre pair (L A , L B ). This low arithmetic complexity algorithm will be utilized as a foundation to search for Legendre pairs with a low time complexity algorithm. The algorithm is presented through the matrix representation of Bluestein's algorithm. Before starting the algorithm, let us define the PSD values. Then, the PSD values of these sequences are defined as: Based on the PSD values of the two sequences, we have the PSD invariant property s.t.
In other words, Legendre pairs are characterized by the constancy of PAF (i.e., the PAF values sum to -2) and the constancy of the PSD (i.e., the PSD values sum to 2ℓ + 2). We illustrate these constancy properties of PAF and PSD with the following example:
Example 3.2. Consider the Legendre pair of order ℓ = 45 (for which m = 22, 2ℓ + 2 = 92)
• the PAF constancy property (with the PAF constant -2) is materialized as follows:
2 -2 s 12 13 14 15 16 17 18 19 20 21 22
• the PSD constancy property (with the PSD constant 2ℓ + 2 = 92) is materialized as follows:
s, PSD(L_A,s), PSD(L_B,s), PSD(L_A,s) + PSD(L_B,s)
1, 18.78398162, 73.21601830, 92 2, 75.60311755, 16.39688232, 92 3, 75.70559038, 16.29440935, 92 4, 67.60652789, 24.39347193, 92 5, 58.80933390, 33.19066596, 92 6, 58.94004336, 33.05995648, 92 7, 16.29533595, 75.70466380, 92 8, 12.42889922, 79.57110066, 92 9, 32.58359208, 59.41640783, 92 10, 8.393304174, 83.60669564, 92 11, 21.17900710, 70.82099280, 92 12, 43.12722512, 48.87277485, 92 13, 72.68711906, 19.31288084, 92 14, 89.20751650, 2.792483344, 92 15, 76, 16, 92 16, 88.23332322, 3.766676806, 92 17, 69.87139244, 22.12860750, 92 18, 59.41640784, 32.58359212, 92 19, 49.23886682, 42.76113297, 92 20, 4.797361782, 87.20263814, 92 21, 6.227140824, 85.77285911, 92 22, 6.864911913, 85.13508805, 92
We note here that the PSD(L A , s) and PSD(L B , s) are the square of the magnitude of the eigenvalues of the matrices A and B generated by the Legendre pair sequences (L A , L B ) in equation ( 19) for s = 1, 2, . . ., ℓ-1 2 .
Remark 3.
where
T n is the symmetric Toeplitz matrix determined via the first column (or row) s.t. t n×1 = t 0 , t 1 , t 2 , . . ., t (ℓ-1) , 0, . . ., 0 t , t j = e -πj 2 i ℓ , and Dℓ = diag(e -πk 2 i ℓ ) ℓ-1 k=0
.
The DFT(L B ) is computed using the same factorization as in (23) with the input sequence L B .
Proof. The symmetric Toeplitz matrix T ℓ is padded with zeros in order to obtain another symmetric Toeplitz matrix T n with a length of n, which is the closest power of 2 to ℓ. The first column (or row) of T n is determined using t n×1 = t 0 , t 1 , t 2 , . . . , t (ℓ-1) , 0, . . . , 0 t . Next, we embed T n into the
, where a symmetric Toeplitz matrix Tn is determined via its first column (or row) s.t. tn×1 = t 0 , 0, . . ., 0, t (ℓ-1) , t (ℓ-2) , . . ., t 2 , t 1 t . Next, we use 2-FFT (as opposed to the well-known 3-FFTs computing a Toeplitz matrix-vector product) to compute
, and F2n is the normalized DFT matrix. Next, we extract
Afterwards, we extract T ℓ from T n using K ℓ×n . Finally, scaling T ℓ by diagonal matrices Dℓ and getting the product of the matrices with each of the two sequences L A and L B of the Legendre pair (L A , L B ) gives the result.
Remark 3.5. 1. The above proposition shows that the matrix representation of Bluestein's FFT algorithm is executed via the symmetric Toeplitz matrices T ℓ [47,44].
We recall here that, computing the Toeplitz matrices by a vector using 2-FFTs [14,37] as opposed to 3-FFTs [9,15] for an even length s.t. n = 2 t had been proposed.
One could also pad the odd length ℓ Legendre pair sequences (L A , L B ) to the closet even length of n = 2 t s.t. n > ℓ, and execute the product of the DFT matrix with each of the two sequences L A and L B of the Legendre pair (L A , L B ) by using the 2-FFT algorithms in [14,37]. Since the purpose is also to search for Legendre pairs of odd length ℓ using a low time complexity-based FFTW in C (as in Section 4), which is the Bluestein's FFT algorithm, we proceed with the matrix factorization of Bluestein's FFT in Proposition 23 in order to reduce arithmetic and time complexities. where n > ℓ.
Proof. Following equation ( 23), the Toeplitz matrix T ℓ is computed using T n followed by 2-FFTs with complexity O(n log(n)). The cost of computing each diagonal matrix Dℓ and D 2n by a vector cost O(ℓ) and O(n), respectively. Thus, the overall cost is O(n log(n)) as n > ℓ.
We recall here that there is also a Winograd FFT algorithm [48] and it minimizes the number of multiplications at the expense of a large number of additions, but this trade-off didn't benefit current processors that need specialized hardware multipliers and was not utilized to implement the FFTW [44].
The PSD shows the influence of noise signals over a spectrum of frequencies. Thus, the PSD serves as a tool for analyzing noise signals and identifying harmonics, guiding the design of filters.
The FFT algorithm is a powerful technique for analyzing the autocorrelation function of a discrete noise signal in order to characterize its PSD. Thus, the FFT-like algorithms are highly efficient in accurately determining the spectral information of dominant noise powers.
In this section, we utilize the proposed FFT-like algorithm with ℓ = 45 for noise filtering with a simple example. We consider a function of time f (t) = cos(ω 1 t)+sin(ω 2 t) with frequencies ω 1 = 40π
and ω 2 = 240π. We then add a large amount of Gaussian white noise to this signal as shown in the first panel of
In this section, we show that the FFTW-based matrix factorization,i.e., Bluestein's algorithm in Section 3 could be utilized to speed up the process of searching Legendre pairs compared to the existing combinatorial searching algorithms. We illustrate these speedups in the particular cases s.t. ℓ = 45 and ℓ = 63. [
Here is a pair of sequences of length 5, that constitutes the 9-compression of the Legendre pair of order 45 presented in Example 3.2 s.t.
[[1, -3, 3, 3, -3], [1, 3, -3, -3, 3]].
Our currently available C code implementations find Legendre pairs of order ℓ = 45 in about 20 minutes, while the FFTW-enabled C codes find Legendre pairs of order ℓ = 45 in less than 10 seconds. Thus, the proposed FFT-like algorithm enables the search for Legendre pairs of order ℓ = 45 with a 99% efficiency, outperforming existing methods.
The factorization ℓ = 63 = 7 • 9 allows one to employ 7-decompression [6] of a pair of sequences of length 9, to construct Legendre pairs of order ℓ = 63. In particular, we employ two copies of the same sequence of length 9 s.t.
This sequence has the following properties: From the above data, one can conclude that the entire search space can be traversed super fast with our novel implementation, while the existing C implementation for Legendre pairs of order ℓ = 63 did not even produce any results i.e. not even one Legendre pair after 24 hours.
In this paper, we provided novel perspectives on Legendre pairs, using concepts from matrix analysis and linear algebra. These perspectives can be used (with great benefits) in current computational schemes to search for Legendre pairs. The structured matrices perspectives on Legendre pairs make use of the Legendre pair matrix equation, to investigate its properties. We introduced a structured matrix approach to obtain a low arithmetic complexity algorithm for computing the product of the DFT matrix with each sequence of the Legendre pair, regardless of their odd lengths. This method enables efficient computation of power spectral density values, and hence, search for Legendre pairs with reduced time complexity. Finally, we showed numerical results based on the C implementation of FFTW in searching for Legendre pairs while attaining a low time complexity algorithm so that the proposed technique excels the conventional combinatorial algorithms in searching Legendre pairs. More specifically, we had shown that the proposed FFT-like algorithm significantly accelerates the search for Legendre pairs of orders 45 and 63, achieving at least 99% improvement in speed compared to conventional algorithms, in the literature.