# Bidiagonal Decompositions
of Vandermonde-Type Matrices of Arbitrary Rank111 This research was partially
supported by Spanish Research Grant PGC2018-096321-B-I00(MCIU/AEI) and by
Gobierno de Aragón (E41-17R) . The authors A. Marco and J. J. Martínez are
members of the Reseach Group asynacs (Ref. ct-ce2019/683) of Universidad de
Alcalá. 222This research was also partially supported by the Woodward Fund for
Applied Mathematics at San José State University. The Woodward Fund is a gift
from the estate of Mrs. Marie Woodward in memory of her son, Henry Teynham
Woodward. He was an alumnus of the Mathematics Department at San José State
University and worked with research groups at NASA Ames.

Jorge Delgado<EMAIL_ADDRESS>Plamen Koev<EMAIL_ADDRESS>Ana Marco
<EMAIL_ADDRESS>José-Javier Martínez<EMAIL_ADDRESS>Juan Manuel
Peña<EMAIL_ADDRESS>Per-Olof Persson<EMAIL_ADDRESS>Steven Spasov
<EMAIL_ADDRESS>
###### Abstract

We present a method to derive new explicit expressions for bidiagonal
decompositions of Vandermonde and related matrices such as the ($q$-, $h$-)
Bernstein-Vandermonde ones, among others. These results generalize the
existing expressions for nonsingular matrices to matrices of arbitrary rank.
For totally nonnegative matrices of the above classes, the new decompositions
can be computed efficiently and to high relative accuracy componentwise in
floating point arithmetic. In turn, matrix computations (e.g., eigenvalue
computation) can also be performed efficiently and to high relative accuracy.

###### keywords:

Vandermonde matrix, totally nonnegative matrix , bidiagonal decomposition ,
eigenvalue

###### MSC:

65F15 , 15A23 , 15B48 , 15B35

††journal: Journal of Computational and Applied Mathematics

[1]organization=Departamento de Matemática Aplicada, Universidad de Zaragoza,
addressline=Edificio Torres Quevedo, city=Zaragoza, postcode=50019,
country=Spain [2]organization=Department of Mathematics, San José State
University, city=San Jose, postcode=95192, country=U.S.A.
[3]organization=Departamento de Física y Matemáticas, Universidad de Alcalá,
addressline=Alcalá de Henares, city=Madrid, postcode=28871, country=Spain
[4]organization=Departamento de Física y Matemáticas, Universidad de Alcalá,
addressline=Alcalá de Henares, city=Madrid, postcode=28871, country=Spain
[5]organization=Departamento de Matemática Aplicada, Universidad de Zaragoza,
addressline=Edificio de Matemáticas, city=Zaragoza, postcode=50019,
country=Spain [6]organization=Department of Mathematics, University of
California, city=Berkeley, postcode=94720, country=U.S.A.
[7]organization=Sofia High School of Mathematics, addressline=61 Iskar St.,
city=Sofia, postcode=1000, country=Bulgaria

##  1 Introduction

A matrix is totally nonnegative (TN) if all of its minors are nonnegative [1,
2, 3]. The bidiagonal decompositions of the TN matrices have become an
important tool in the study of these matrices [4, 5] and for performing matrix
computations with them accurately and efficiently [6, 7, 8]. The bidiagonal
decompositions of many classical TN matrices such as Vandermonde and many
related matrices are well known, but contain singularities when some of the
nodes coincide. For example, a $3\times 3$ Vandermonde matrix with nodes
$x,y,z$ is decomposed as [6]:

$\displaystyle\left[\begin{array}[]{ccc}1&x&x^{2}\\\ 1&y&y^{2}\\\
1&z&z^{2}\end{array}\right]$ $\displaystyle=\left[\begin{array}[]{ccc}1\\\
&1\\\ &1&1\end{array}\right]\left[\begin{array}[]{ccc}1\\\ 1&1\\\
&\frac{z-y}{y-x}&1\end{array}\right]\left[\begin{array}[]{ccc}1\\\ &y-x\\\
&&(z-x)(z-y)\end{array}\right]$ (13) $\displaystyle\phantom{=}\mbox{
}\;\times\left[\begin{array}[]{ccc}1&x\\\ &1&y\\\
&&1\end{array}\right]\left[\begin{array}[]{ccc}1\\\ &1&x\\\
&&1\end{array}\right].$ (20)

This decomposition is undefined when $x=y$, which is unfortunate, since the
Vandermonde matrix is very well defined for any values of the nodes.

By relaxing the requirement that the bidiagonal factors have ones on the main
diagonal, we show how to rearrange the factors in (20), so that the new
decomposition contains no singularities and is valid for any values of the
nodes and is thus valid for a Vandermonde matrix of arbitrary rank. For
example, the $3\times 3$ Vandermonde from (20) can be decomposed as

$\displaystyle\left[\begin{array}[]{ccc}1&x&x^{2}\\\ 1&y&y^{2}\\\
1&z&z^{2}\end{array}\right]$ $\displaystyle=\left[\begin{array}[]{ccc}1\\\
&1\\\ &1&z-y\end{array}\right]\left[\begin{array}[]{ccc}1\\\ 1&y-x\\\
&1&z-x\end{array}\right]\left[\begin{array}[]{ccc}1\\\ &1\\\
&&1\end{array}\right]$ (33) $\displaystyle\phantom{=}\mbox{
}\;\times\left[\begin{array}[]{ccc}1&x\\\ &1&y\\\
&&1\end{array}\right]\left[\begin{array}[]{ccc}1\\\ &1&x\\\
&&1\end{array}\right],$ (40)

which is valid for any values of the nodes $x,y,z$.

Many classes of other TN matrices share the exact same singularities in their
bidiagonal decompositions, e.g., the Bernstein-Vandermonde, their $q$-, $h$-,
and rational generalizations, Lupaş, Said-Ball matrices, etc. [6, 9, 10, 11,
12, 13, 14, 15, 16]. We call these matrices Vandermonde-type below (see
section 3).

The method that allowed us to obtain the decomposition (40) from (20) applies
to all Vandermonde-type matrices (section 5) and is the main contribution of
this paper. The starting point for our method is the existing bidiagonal
decompositions of the Vandermonde-type matrices, which are valid only when
these matrices are both TN and nonsingular. Our method takes these
decompositions as a starting point and produces new bidiagonal decompositions
valid for matrices of arbitrary rank, regardless of whether they are TN or not
– see Corollary 4.1.

While the computational complexity of the transformed bidiagonal decomposition
of an $n\times n$ matrix remains $O(n^{2})$, the new expressions are simpler.

In terms of accuracy, just like their non-singular counterparts, the new
decompositions remain insusceptible to subtractive cancellation, and thus all
of the entries of these decompositions can be computed to high relative
accuracy when the matrix is TN. By “high relative accuracy” we mean that for
each entry its sign and most of its leading significant digits are computed
correctly (see section 8).

The Vandermonde-type matrices have many applications that are well referenced
in the papers we cite above that deal with the nonsingular case. For example,
the Lupaş matrices (section 10.1) have direct applications in CAGD [11, 17].
For TN Vandermonde-type matrices, the new results in this paper allow for
matrix computations with them to be performed to high relative accuracy very
efficiently (in $O(n^{3})$ time) using the methods of [8] now also when these
matrices are singular – see section 9 for a numerical example.

The efficiency and high relative accuracy is particularly relevant, for
example, in eigenvalue computations since the corresponding matrices are
unsymmetric. The error bounds for the eigenvalues computed by the conventional
algorithms (such as the ones in LAPACK [18]) [19] imply that none of the
eigenvalues are guaranteed to be accurate, although the largest ones typically
are – see the example in section 9. In contrast, the results of this paper now
allow for all eigenvalues to be efficiently computed to high relative accuracy
and, in particular, the zero eigenvalues are computed exactly.

The paper is organized as follows. In section 2 we review the bidiagonal
decompositions of nonsingular TN matrices. In section 3 we present the formal
definition of Vandermonde-type matrix and, in section 4, we show how to remove
the singularities in its bidiagonal decomposition. We demonstrate how our
method applies to the particular class of $q$-Bernstein-Vandermonde matrices
in section 5. Our method is also directly applicable to derivative matrices,
such as generalized Vandermonde matrices, Laguerre matrices, etc., which are
submatrices or products of Vandermonde matrices and other nonsingular TN
matrices (section 6). In section 7, we present a method to make the bottom
right-hand corner entry of all bidiagonal factors in the corresponding
bidiagonal decompositions equal to 1 – this is a requirement for the methods
of [8] to work. We discuss accuracy issues in section 8 and present numerical
experiments in section 9. The explicit formulas for the decompositions of
several Vandermonde-type matrices are in the Appendix.

##  2 Bidiagonal decompositions of TN matrices

Our focus in this paper is on the class of TN matrices, but the formulas and
methods we present are valid without the requirement of total nonnegativity
(Corollary 4.1 below). The bidiagonal decompositions of the TN matrices serve
as a major tool in their study and computations with them. We review those
here [2, 4].

A nonsingular $n\times n$ TN matrix $A$ can be uniquely factored as

$A=L^{(1)}L^{(2)}\cdots L^{(n-1)}DU^{(n-1)}U^{(n-2)}\cdots U^{(1)},$ (41)

where $L^{(i)}$ are $n\times n$ nonnegative and unit lower bidiagonal, $D$ is
$n\times n$ nonnegative and diagonal, and $U^{(i)}$ are $n\times n$
nonnegative and unit upper bidiagonal. For the nontrivial entries of the
factors $L^{(k)}$ and $U^{(k)}$, respectively, we have
$L_{i+1,i}^{(k)}=U_{i,i+1}^{(k)}=0$ for $i<n-k$.

We call the above decomposition (41), the ordinary bidiagonal decomposition of
$A$ to distinguish it from what we will define below as the singularity-free
bidiagonal decomposition of a TN matrix.

The decomposition (41) occurs naturally in the process of complete Neville
elimination when adjacent rows and columns are used for elimination [4, 5].

There are exactly $n^{2}$ nontrivial entries which parameterize the ordinary
bidiagonal decomposition (41). Following [6, sec. 4], those nontrivial entries
can be conveniently stored in an $n\times n$ matrix $M$, denoted
$M={\mathcal{B}D}(A)$, where:

$\displaystyle m_{ij}$ $\displaystyle=L_{i,i-1}^{(n-i+j)},\quad i>j,$
$\displaystyle m_{ij}$ $\displaystyle=U_{j-1,j}^{(n-j+i)},\quad i<j,$ (42)
$\displaystyle m_{ii}$ $\displaystyle=D_{ii}.$

Namely, for $k=1,2,\ldots,n-1$,

$L^{(k)}=\left[\begin{array}[]{ccccccc}1\\\ &\ddots\\\ &m_{n-k+1,1}&1\\\
&&m_{n-k+2,2}&1\\\ &&&\ddots&\ddots\\\ &&&&m_{nk}&1\end{array}\right],$
$U^{(k)}=\left[\begin{array}[]{ccccccc}1\\\ &\ddots\\\ &&1&m_{1,n-k+1}\\\
&&&1&m_{2,n-k+2}\\\ &&&&\ddots&\ddots\\\ &&&&&1&m_{nk}\\\
&&&&&&1\end{array}\right],$

and

$D=\mbox{diag}(m_{11},m_{22},\ldots,m_{nn}).$

For $i\neq j$, the $m_{ij}$ are the multipliers of the complete Neville
elimination with which the $(i,j)$ entry of $A$ is eliminated and $m_{ii}$ are
the diagonal entries of $D$.

For example, if $A$ is the $3\times 3$ Vandermonde matrix from (20), then the
matrix $M={\mathcal{B}D}(A)$ is:

$M=\left[\begin{array}[]{ccc}1&x&x\\\ 1&y-x&y\\\
1&\frac{z-y}{y-x}&(z-x)(z-y)\end{array}\right].$ (43)

In this paper, we assume that explicit expressions for the entries of the
ordinary bidiagonal decomposition (41) are given. For the matrices we consider
in this paper, those expressions contain singularities when some of the nodes
coincide and we work to remove those singularities.

To that end, we drop the requirement that the matrices $L^{(i)}$ and $U^{(i)}$
have ones on the main diagonal. We define a new bidiagonal decomposition,
which we call singularity-free bidiagonal decomposition and denote it by
$\mathcal{SBD}(A)$, as

$A=L_{1}L_{2}\cdots L_{n-1}DU_{n-1}U_{n-2}\cdots U_{1}.$ (44)

The matrix $D$ is nonnegative and diagonal. The factors $L_{i}$ and $U_{i}$
are nonnegative lower and upper bidiagonal, and have the same nonzero patterns
as $L^{(i)}$ and $U^{(i)}$, respectively, $i=1,2,\ldots,n-1$. Namely,
$(L_{k})_{i+1,i}=(U_{k})_{i,i+1}=0$ for $i<n-k$ (see theorem 2.1 in [6]).

Following [8], the nontrivial entries of $\mathcal{SBD}(A)$ are stored in two
matrices: $B$, which is $n\times n$ and $C$, which is $(n+1)\times(n+1)$:

$[B,C]=\mathcal{SBD}(A).$

As with $\mathcal{BD}(A)$, the matrix $B$ stores the nontrivial offdiagonal
entries of $L_{i}$ and $U_{i}$ as well as the diagonal entries of $D$, exactly
as in (42):

$\displaystyle b_{ij}$ $\displaystyle=(L_{n-i+j})_{i,i-1},i>j,$ $\displaystyle
b_{ij}$ $\displaystyle=(U_{n-j+i})_{j-1,j},i<j,$ $\displaystyle b_{ii}$
$\displaystyle=D_{ii}.$

The matrix $C$ stores the diagonal entries of $L_{i}$ and $U_{i}$ as

$c_{ij}=\left\\{\begin{array}[]{ll}(L_{n-i+j})_{i-1,i-1},&i>j,\\\
(U_{n-j+i})_{j-1,j-1},&i<j.\end{array}\right.$

In this arrangement, $c_{ij},\,i>j$, is the diagonal entry in $L_{n-i+j}$
immediately above $b_{ij}$ and similarly for $i<j$ and $U_{n-j+i}$. The
entries $c_{ii},i=1,2,\ldots,n+1$ as well as the entries $c_{1,n+1}$ and
$c_{n+1,1}$ are unused. This is the same construction as the one given in
formula (9) in [8], except that we now allow for the $(n,n)$ entry in $L_{i}$
and $U_{i}$ to be any nonnegative number and not necessarily equal 1. This is
the reason we need an $(n+1)\times(n+1)$ matrix to host the entries $c_{ij}$:
the nontrivial diagonal entries of $L_{1},L_{2},\ldots,L_{n-1}$ are of lengths
$2,3,\ldots,n$, and similarly for the $U_{i}$’s.

As we explain in section 7, we can always make the $(n,n)$ entry in the
$L_{i}$’s and $U_{i}$’s equal to 1, so that the algorithms of [8] be used, but
the formulas are more elegant without that restriction, which can be imposed
after in software.

Thus, in this paper we start with the known ordinary bidiagonal decomposition
$M=\mathcal{BD}(A)$ of a nonsingular matrix $A$, (41), and produce the
singularity-free bidiagonal decomposition $[B,C]=\mathcal{SBD}(A)$, (44),
which is defined for a matrix $A$ of arbitrary rank.

The new singularity-free bidiagonal decomposition (44) is not unique, but this
is inconsequential – for the purposes of accurate computations, any accurate
decomposition of a TN matrix as a product of nonnegative bidiagonals is an
equally good input [6].

##  3 Vandermonde-type matrices

The ordinary Vandermonde matrices are the starting point of our investigation
as the approach extends analogously to all other classes of matrices in this
paper.

An $n\times n$ Vandermonde matrix $A$ with nodes $x_{1},x_{2},\ldots,x_{n}$ is
defined as

$A=\left[x_{i}^{j-1}\right]_{i,j=1}^{n}.$ (45)

When the nodes are distinct, it has an ordinary bidiagonal decomposition (41),
$V={\cal BD}(A)$, such that [6]

$\begin{array}[]{lcll}v_{ii}&=&\displaystyle\prod_{k=1}^{i-1}(x_{i}-x_{k}),&\quad
1\leq i\leq n,\\\
v_{ij}&=&\displaystyle\prod_{k=i-j}^{i-2}\frac{x_{i}-x_{k+1}}{x_{i-1}-x_{k}},&\quad
1\leq j<i\leq n,\\\ v_{ij}&=&x_{i},&\quad 1\leq i<j\leq n.\end{array}$ (46)

In the following definition, we take a more general stance, where the lower
bidiagonal factors and the diagonal have additional factors, which contain no
singularities.

######  Definition 1.

An $n\times n$ matrix $A$ is said to be of Vandermonde-type with nodes
$x_{1},x_{2},\ldots,x_{n}$ on an (open or closed) interval $D$, if it
satisfies all of the conditions below:

  1. 1.

$a_{ij}=f_{j}(x_{i})$, where $f_{j}(x)$ is a rational function of $x$ for
$i,j=1,2,\ldots,n$;

  2. 2.

$A$ is TN when $x_{1}<x_{2}<\cdots<x_{n}$ and all $x_{i}\in
D,\;i=1,2,\ldots,n$;

  3. 3.

$A$ has an ordinary bidiagonal decomposition, $M={\cal BD}(A)$, such that

$m_{ij}=s_{ij}v_{ij},$ (47)

for $1\leq j\leq i\leq n$, where $v_{ij}$ are defined as in (46), and the
entries $s_{ij}$ for $1\leq j\leq i\leq n$ and the entries $m_{ji}$ for $1\leq
j<i\leq n$ are rational functions of $x_{1},x_{2},\ldots,x_{n}$, with no
singularities when $x_{i}\in D,\,i=1,2,\ldots,n$.

A Vandermonde-type matrix will thus be fully specified by the entries $s_{ij}$
for $1\leq j\leq i\leq n$ and the entries $m_{ji}$ for $1\leq j<i\leq n$.

From the above definition, the ordinary Vandermonde matrix (45) is a
Vandermonde-type matrix on $(0,\infty)$ with $s_{ij}=1$ for $1\leq j\leq i\leq
n$ and $m_{ji}=v_{ji}=x_{j}$ for $1\leq j<i\leq n$.

In the literature, the Vandermonde-type matrices sometimes have $n+1$ nodes,
with the indexing starting at 0, e.g., $x_{0},x_{1},\ldots,x_{n}$. For
consistency, in this paper, all Vandermonde-type matrices are $n\times n$ and
their nodes are $x_{1},x_{2},\ldots,x_{n}$.

##  4 Singularity-free bidiagonal decompositions of Vandermonde-type matrices

In this section, we show how to remove the singularities in the ordinary
bidiagonal decomposition (41) for Vandemonde-type matrices and obtain a
singularity-free bidiagonal decomposition (44).

If $M={\mathcal{B}D}(A)$ is the ordinary bidiagonal decomposition (41) of a
Vandermonde-type matrix, from (47) we have
$m_{ii}=s_{ii}v_{ii},\,i=1,2,\ldots,n$, and thus the diagonal factor $D$ can
be factored as $D=D_{n}\cdot D^{\prime}$, so that

$A=L^{(1)}L^{(2)}\cdots L^{(n-1)}\cdot D_{n}\cdot D^{\prime}\cdot
U^{(n-1)}\cdots U^{(1)},$

where

$D_{n}={\rm diag}(v_{11},\ldots,v_{nn})$ (48)

and $D^{\prime}={\rm diag}(s_{11},\ldots,s_{nn})$.

Since there are no singularities in the product $D^{\prime}\cdot
U^{(n-1)}\cdots U^{(1)}$, it suffices to remove the singularities in the
remaining factors, i.e., in the product

$L^{(1)}L^{(2)}\cdots L^{(n-1)}D_{n}.$ (49)

######  Theorem 4.1.

Let $L^{(i)},\,i=1,2,\ldots,n-1$ be the lower bidiagonal factors in the
ordinary bidiagonal decomposition (41) of a Vandermonde-type matrix, whose
nontrivial entries are defined as in (47) and $D_{n}$ be defined as in (48).
Then

$L_{1}L_{2}\cdots L_{n-1}=L^{(1)}L^{(2)}\cdots L^{(n-1)}D_{n},$

where

$L_{i}\\!=\\!\left[\begin{array}[]{ccccccc}1\\\ &\ddots\\\
&&\\!\\!\\!\\!\\!\\!\\!\\!1\\\
&&\\!\\!\\!\\!\\!\\!\\!\\!s_{n-i+1,1}&x_{n-i+1}-x_{n-i}\\\
&&&s_{n-i+2,2}&x_{n-i+2}-x_{n-i}\\\ &&&&\\!\\!\\!\\!\\!\\!\\!\ddots&\ddots\\\
&&&&&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!s_{ni}&\\!\\!\\!\\!\\!\\!\\!\\!x_{n}-x_{n-i}\end{array}\right]\\!.$
(50)

###### Proof.

It suffices to observe just the initial step, i.e., that

$D_{n-1}\cdot L_{n-1}=L^{(n-1)}\cdot D_{n},$ (51)

where $D_{n-1}$ is a diagonal matrix, such that $(D_{n-1})_{ii}=1$ for
$i=1,2,$ and

$(D_{n-1})_{ii}=\prod_{k=2}^{i-1}(x_{i}-x_{k}),\quad i=3,4,\ldots,n.$

In other words, the $(n-1)\times(n-1)$ bottom right principal submatrix of
$D_{n-1}$ is the diagonal factor in the ordinary bidiagonal decomposition of
an $(n-1)\times(n-1)$ ordinary Vandermonde matrix with nodes
$x_{2},x_{3},\ldots,x_{n}$.

Since both $D_{n}$ and $D_{n-1}$ are diagonal and both $L^{(n-1)}$ and
$L_{n-1}$ are lower bidiagonal, we have bidiagonals on each side of (51). To
establish that those are equal, it suffices to show that the corresponding
diagonal and offdiagonal entries are the same. Since $L^{(n-1)}$ is unit lower
bidiagonal, the diagonal entry in the product on the right, $L^{(n-1)}D_{n}$
is $(D_{n})_{ii}$. Thus, for $i>1$,

$(D_{n})_{ii}=\prod_{k=1}^{i-1}(x_{i}-x_{k})=\left(\prod_{k=2}^{i-1}(x_{i}-x_{k})\right)\cdot(x_{i}-x_{1})=(D_{n-1})_{ii}(L_{n-1})_{ii},$

since $(L_{n-1})_{ii}=x_{i}-x_{1}$ by (50). Also, $(D_{n})_{11}=1$, so the
diagonal entries on both sides of (51) are equal.

Since $(L_{n-1})_{i+1,i}=s_{i+1,i}$, the offdiagonal entry in position
$(i+1,i)$ on the left side of (51) is $s_{i+1,i}(D_{n-1})_{i+1,i+1}$. Also,
from (42),

$L_{i+1,i}^{(n-1)}=m_{i+1,i}=s_{i+1,i}v_{i+1,i}=s_{i+1,i}\prod_{j=1}^{i-1}\frac{x_{i+1}-x_{j+1}}{x_{i}-x_{j}}$

and thus

$\displaystyle s_{i+1,i}(D_{n-1})_{i+1,i+1}$
$\displaystyle=s_{i+1,i}\prod_{j=2}^{i}(x_{i+1}-x_{j})$
$\displaystyle=s_{i+1,i}\prod_{j=1}^{i-1}(x_{i+1}-x_{j+1})$
$\displaystyle=s_{i+1,i}\prod_{j=1}^{i-1}\frac{x_{i+1}-x_{j+1}}{x_{i}-x_{j}}\prod_{j=1}^{i-1}(x_{i}-x_{j})$
$\displaystyle=s_{i+1,i}v_{i+1,i}(D_{n})_{ii}$
$\displaystyle=(L^{(n-1)})_{i+1,i}\,(D_{n})_{ii},$

which is the $(i+1,i)$ entry on the right hand side of (51).

Therefore the offdiagonal entries on each side of (51) are also equal and (51)
is fully established.

Completely analogously, we establish that $D_{k-1}L_{k-1}=L^{(k-1)}D_{k}$ for
$k=n-1,n-2,\ldots,2,$ and for the product (49) we have

$\displaystyle L^{(1)}L^{(2)}\cdots\underline{L^{(n-1)}D_{n}}$
$\displaystyle=L^{(1)}L^{(2)}\cdots\underline{L^{(n-2)}D_{n-1}}L_{n-1}$
$\displaystyle=L^{(1)}L^{(2)}\cdots\underline{L^{(n-3)}D_{n-2}}L_{n-2}L_{n-1}$
$\displaystyle=\cdots$ $\displaystyle=\underline{L^{(1)}D_{2}}L_{2}\cdots
L_{n-1}$ $\displaystyle=\underline{D_{1}L_{1}}L_{2}\cdots L_{n-1},\quad\mbox{
and since $D_{1}=I$},$ $\displaystyle=L_{1}L_{2}\cdots L_{n-1}.$

The factors that change on each step are underlined.

In other words, we “walk” the matrix $D_{n}$ through the product (49) right-
to-left, canceling the denominators in each $L^{(k)}$ and factoring a slightly
different diagonal matrix, $D_{k-1}$, out the left until we end up with just
$D_{1}=I$, which we drop. ∎

The ordinary bidiagonal decompositions exist for nonsingular TN matrices,
however very simple singular TN matrices such as
$\big{[}\genfrac{}{}{0.0pt}{}{0\;\;0}{1\;\;0}\big{]}$ (let alone non-TN ones)
don’t have ordinary bidiagonal decompositions. Thus the fact that all
Vandermonde-type matrices have singularity-free bidiagonal decompositions
regardless of their rank or total nonnegativity is remarkable, which we prove
below.

######  Corollary 4.1.

The singularity-free bidiagonal decomposition (44) of a Vandermonde-type
matrix $A$ is valid for any complex values of the nodes
$x_{1},x_{2},\ldots,x_{n}$ for which the matrices on both sides of (44) are
defined.

###### Proof.

Since it is derived from the ordinary bidiagonal decomposition (41), the
singularity-free bidiagonal decomposition (44) is valid when $A$ is
nonsingular and TN. This remains the case when all the nodes of $A$ are
strictly increasing inside an open interval $(a,b)$ contained in the interval
$D$ of total nonnegativity of $A$.

The $(i,j)$th entries on each side of (44) are rational functions of $x_{i}$
for all $i,j=1,2,\ldots,n$. Indeed, on the left hand side,
$a_{ij}=f_{j}(x_{i})$, is a rational function of $x_{i}$ by definition 1. On
the right hand side, the $(i,j)$th entry is a polynomial in the entries of the
factors of (44). These entries are obtained via the construction in section 4
from the entries of the ordinary bidiagonal decomposition (41), which are
rational functions of the nodes $x_{1},x_{2},\ldots,x_{n}$ as either quotients
of minors of $A$ or products of quotients of minors of $A$ [6, Prop. (3.1)].

Since all rational functions are meromorphic on $\mathbb{C}$, and the equality
between the $(i,j)$th entries on each side of (44) holds on an open interval
in $\mathbb{R}$ containing $x_{i}$: $x_{i-1}<x_{i}<x_{i+1}$, for $1<i<n$ and
$a<x_{1}<x_{2}$, $x_{n-1}<x_{n}<b$ for $i=1,i=n$, respectively, the Identity
Theorem [20, Thm. 3.2.6], implies that this equality holds for any complex
values of the nodes where these entries are defined. ∎

As a direct application to the ordinary Vandemonde matrices for example, as is
evident from the formulas in section 5.1, the above Corollary immediately
implies that their singularity-free bidiagonal decomposition is valid for any
complex nodes .

##  5 Our method and explicit formulas for Vandemonde-type matrices

Our method for producing a singularity-free bidiagonal decomposition for an
$n\times n$ Vandermonde-type matrix $A$ of arbitrary rank works as follows:

  1. 1.

We start with the ordinary bidiagonal decomposition $M=\mathcal{BD}(A)$. Since
$A$ is Vandermonde-type, the entries $m_{ij}$ factor as:

$m_{ij}=s_{ij}v_{ij}$

for $1\leq j\leq i\leq n$, where $v_{ij}$ are defined in (46) and $s_{ij}$
contain no singularities.

  2. 2.

The singularity-free bidiagonal decomposition $[B,C]=\mathcal{SBD}(A)$, is
then

$\displaystyle b_{ij}$ $\displaystyle=s_{ij},$ $\displaystyle\quad 1\leq j\leq
i\leq n,$ (52) $\displaystyle b_{ij}$ $\displaystyle=m_{ij},$
$\displaystyle\quad 1\leq i<j\leq n,$ (53) $\displaystyle c_{ij}$
$\displaystyle=x_{i-1}-x_{i-j},$ $\displaystyle\quad 2\leq j<i\leq n+1,$ (54)
$\displaystyle c_{ij}$ $\displaystyle=1,$ $\displaystyle\quad{j>i}.$ (55)

In other words, to obtain the lower bidiagonal factors $L_{i}$ in (44) from
the corresponding factors $L^{(i)}$ in (41), we remove the factors $v_{ij}$
from subdiagonal elements in $L^{(i)}$ and set the diagonal entries in
$L^{(i)}$ as in (54).

To obtain the diagonal factor in (44) we remove the factors $v_{ii}$ from the
corresponding diagonal entry $m_{ii}$ in the diagonal factor in (41),
$i=1,2,\ldots,n$.

The upper bidiagonal factors $U^{(i)}$ remain intact,
$U_{i}=U^{(i)},\,i=1,2,\ldots,n-1,$ as per (53) and (55).

To fully specify the entire decomposition, $[B,C]=\mathcal{SBD}(A)$, it
suffices, therefore, to specify the entries $s_{ij},i\geq j$, and
$m_{ij},i<j$, The rest of the entries in $B$ and $C$ are specified by (53) and
(55).

###  5.1 Ordinary Vandermonde matrices

The ordinary bidiagonal decomposition of an ordinary Vandermonde matrix is
given in (46).

Their singularity-free bidiagonal decomposition is such that

$b_{ij}=1$

for $i\geq j$ and

$b_{ij}=x_{j}$

for $i<j$ and $c_{ij}=x_{i-1}-x_{i-j}$ for $2\leq j<i\leq n+1$ and $c_{ij}=1$
otherwise. In other words, the factors in (44) are

$L_{i}\\!=\\!\left[\begin{array}[]{ccccccc}1\\\ &\ddots\\\ &&1\\\
&&1&x_{n-i+1}-x_{n-i}\\\ &&&1&x_{n-i+2}-x_{n-i}\\\
&&&&\\!\\!\\!\\!\\!\\!\\!\ddots&\ddots\\\
&&&&&\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!1&\\!\\!\\!\\!\\!\\!\\!\\!x_{n}-x_{n-i}\end{array}\right],$

$D=I$ and

$U_{i}=\left[\begin{array}[]{ccccccc}1\\\ &\ddots\\\ &&1&x_{n-i}\\\
&&&\ddots&\ddots\\\ &&&&1&x_{n-1}\\\ &&&&&1\\\ \end{array}\right].$

The singularity-free bidiagonal decomposition of a an ordinary Vandermonde
matrix can be computed to high relative accuracy in $O(n^{2})$ time using the
routine `STNBDVandermonde` in our package STNTool [21].

For example, if $A$ is the $3\times 3$ Vandermonde matrix in our example (40)
from the Introduction, then its singularity-free bidiagonal decomposition
$[B,C]={\mathcal{S}BD}(A)$ is such that

$B=\left[\begin{array}[]{ccc}1&x&x\\\ 1&1&y\\\
1&1&1\end{array}\right]\quad\mbox{and}\quad
C=\left[\begin{array}[]{cccc}1&1&1&1\\\ 1&1&1&1\\\ 1&y-x&1&1\\\
1&z-y&z-x&1\end{array}\right].$

The same technique applies to all Vandemonde-type matrices. We demonstrate the
derivation of one non-trivial case, the $q$-Bernstein Vandemonde matrices. In
the Appendix (section 10), we list the explicit formulas for several other
classes of Vandermonde-type matrices.

###  5.2 $q$-Bernstein-Vandemonde matrices

The $q$-Bernstein-Vandermonde matrices are expressed in terms of $q$-integers
and $q$-binomial coefficients: given $q>0$ and any nonnegative integer $r$, a
$q$-integer $[r]$ is defined as

$[r]=\left\\{\begin{array}[]{ll}(1-q^{r})/(1-q),&q\neq 1,\\\
r,&q=1.\end{array}\right.$ (56)

For a nonnegative integer $r$, a $q$-factorial is defined as

$[r]!=\left\\{\begin{array}[]{ll}[r][r-1]\cdots[1],&r\neq 1,\\\
1,&r=0.\end{array}\right.$

The $q$-binomial coefficient is defined as

$\left[\genfrac{}{}{0.0pt}{0}{n}{r}\right]=\frac{[n][n-1]\cdots[n-r+1]}{[r]!}$

for integers $n\geq r\geq 0$ and as zero otherwise.

The $q$-Bernstein-Vandermonde matrices generalize the Bernstein-Vandermonde
matrices [12, 14], which are the $q=1$ case. The nonsingular case was
developed in [10].

The $q$-Bernstein polynomials of degree $n$ for $0<q\leq 1$ are defined in
[22] as

$b_{i,q}^{n}(x)=\left[\genfrac{}{}{0.0pt}{0}{n}{i}\right]x^{i}\prod_{s=0}^{n-i-1}(1-q^{s}x),\quad
x\in[0,1],\quad i=0,1,\ldots,n.$ (57)

An $n\times n$ $q$-Bernstein-Vandermonde matrix $A$ with nodes
$x_{1},x_{2},\ldots,x_{n}$ is defined as

$A=\left[b^{n-1}_{j-1,q}(x_{i})\right]_{i,j=1}^{n}$ (58)

and is TN when

$0\leq x_{1}\leq x_{2}\leq\cdots\leq x_{n}<1.$ (59)

When the nodes $x_{i},\;i=1,2,\ldots,n$, are also distinct, it has an ordinary
bidiagonal decomposition $M=\mathcal{BD}(A)$ with parameters [10]

$\displaystyle m_{ij}$
$\displaystyle=\frac{1-q^{n-j}x_{i-j}}{1-q^{n-j}x_{i-1}}\prod_{s=0}^{n-1-j}\frac{1-q^{s}x_{i}}{1-q^{s}x_{i-1}}\underbrace{\prod_{k=i-j}^{i-2}\frac{x_{i}-x_{k+1}}{x_{i-1}-x_{k}}}_{v_{ij}},$
$\displaystyle m_{ji}$
$\displaystyle=\frac{[n-i+1]}{[i-1]}\cdot\frac{x_{j}}{1-q^{n-i}x_{j}}\prod_{k=1}^{j-1}\frac{1-q^{n-i+1}x_{k}}{1-q^{n-i}x_{k}},$
(60)

for $1\leq j<i\leq n$ and

$\displaystyle m_{ii}$
$\displaystyle=\left[\genfrac{}{}{0.0pt}{0}{n-1}{i-1}\right]\frac{\prod_{s=0}^{n-1-i}(1-q^{s}x_{i})}{\prod_{k=1}^{i-1}(1-q^{n-i}x_{k})}\underbrace{\prod_{k=1}^{i-1}(x_{i}-x_{k})}_{v_{ii}}$

for $1\leq i\leq n$.

With the entries of the ordinary bidiagonal decomposition of the ordinary
Vandermonde matrix $v_{ij}$ defined as in (46), we have that
$m_{ij}=s_{ij}v_{ij}$ for $1\leq j\leq i\leq n$, where

$\displaystyle s_{ij}$
$\displaystyle=\frac{1-q^{n-j}x_{i-j}}{1-q^{n-j}x_{i-1}}\prod_{s=0}^{n-1-j}\frac{1-q^{s}x_{i}}{1-q^{s}x_{i-1}},$
(61) $\displaystyle s_{ii}$
$\displaystyle=\left[\genfrac{}{}{0.0pt}{0}{n-1}{i-1}\right]\frac{\prod_{s=0}^{n-1-i}(1-q^{s}x_{i})}{\prod_{k=1}^{i-1}(1-q^{n-i}x_{k})}.$
(62)

The $q$-Bernstein-Vandermonde matrix is therefore a Vandermonde-type matrix on
$[0,1)$.

As described at the beginning of this section, we obtain the singularity-free
bidiagonal decomposition of $A$ by removing the factors $v_{ij}$ from $m_{ij}$
for $i\geq j$ and setting $c_{ij}=x_{i-1}-x_{i-j}$ for $2\leq j<i\leq n+1$.

The decomposition $[B,C]=\mathcal{SBD}(A)$ is thus fully defined as in
(52)–(55) with the $s_{ij},\,i\leq j$, defined in (61) and (62) and
$m_{ji},\,i>j,$ defined in (60).

The singularity-free bidiagonal decomposition of a $q$-Bernstein-Vandermonde
matrix can be computed to high relative accuracy in $O(n^{2})$ time using the
routine `STNBDqBernsteinVandermonde` in our package STNTool [21].

##  6 Derivative Vandermonde matrices

Our results apply to other classes of matrices which are not directly
Vandermonde-type. These include

  * 1.

submatrices of Vandermonde-type matrices,

  * 2.

products of Vandermonde-type matrices and other matrices, whose ordinary
bidiagonal decompositions are known and contain no singularities.

The former include the generalized Vandermonde matrices. The latter include
the Laguerre matrices, the Bessel matrices, and Wronskian matrices of a basis
of exponential polynomials. We address these separately.

The generalized Vandermonde matrices

$G_{\lambda}(x_{1},\ldots,x_{n})=\left[x_{i}^{j-1+\lambda_{n-j+1}}\right]_{i,j=1}^{n}$

are defined for integer partitions
$\lambda=(\lambda_{1},\lambda_{2},\ldots,\lambda_{n})$, where $\lambda_{i}$
are integers such that
$\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}\geq 0$. They are
submatrices of ordinary Vandermonde matrices (obtained by removing appropriate
columns). Their singularity-free bidiagonal decompositions can thus be
obtained by starting with the singularity-free bidiagonal decomposition of the
ordinary Vandermonde matrix, derived in section 5.1, and removing the
appropriate columns employing the methods of [8].

The singularity-free bidiagonal decomposition of a generalized Vandermonde
matrix can be computed to high relative accuracy in $O(\lambda_{1}n^{2})$ time
using the `STNBDGeneralizedVandermonde` routine in our package STNTool [21].

Several classes of TN matrices are products of a TN ordinary Vandermonde
matrix (call it $V$) and a second TN matrix (call it $A$) whose ordinary
bidiagonal decomposition is known and contains no singularities. These
include:

  * 1.

Laguerre [23],

  * 2.

Bessel matrices [24], and

  * 3.

Wronskian matrices of a basis of exponential polynomials [25].

The singularity-free bidiagonal decomposition of the product $VA$ (or $AV$)
can be obtained using the method for computing the singularity-free bidiagonal
decomposition of a product of TN matrices from [8, sec. 6] (implemented in the
routine `STNProduct` in our package STNTool [21]) with the singularity-free
bidiagonal decomposition of $V$ from section 3 and the ordinary bidiagonal
decomposition of $A$.

##  7 The $(n,n)$ entry of the bidiagonal factors must equal 1 for
computations

The algorithms of [8] require that the $(n,n)$ entry of every bidiagonal
factor in the singularity-free bidiagonal decomposition of a TN matrix be
equal to 1. These algorithms rely heavily on this fact and it turns out that
the assumption that this is always the case can be made without any loss of
generality.

To that end, in this section, we present a method to “fix” any singularity-
free bidiagonal decomposition (44) so that the bottom right-hand corner
entries of all bidiagonal factors, $L_{1},L_{2},\ldots,L_{n-1}$ and
$U_{1},U_{2},\ldots,U_{n-1}$ equal to 1.

We do so by “moving” the non-unit $(n,n)$ entries of each lower bidiagonal
factor into $D$ and adjusting the $(n,n-1)$ entries accordingly:

For any bidiagonal matrix $L_{i}$ and diagonal matrix $D_{i}={\rm
diag}(1,1,\ldots,1,x_{i}),$ $i=1,2,\ldots,n-1,$ we can make the bottom right
hand corner entry of the product $D_{i}L_{i}$ equal to one by factoring a new
diagonal factor $D_{i+1}$ out of the right. Namely,

$D_{i}L_{i}=\overline{L}_{i}D_{i+1},$ (63)

where $\overline{L}_{i}$ equals $L_{i}$, except for
$(\overline{L}_{i})_{nn}=1,\,(\overline{L}_{i})_{n,n-1}=L_{n,n-1}x_{i}$. Then
$D_{i+1}={\rm diag}(1,1,\ldots,x_{i}L_{nn})$.

This allows us, by starting with $D_{1}=I$, to move the bottom right hand
corner entries of the lower bidiagonal factors of $\mathcal{SBD}(A)$ into the
diagonal factor, $D$:

$\displaystyle A$ $\displaystyle=L_{1}L_{2}\cdots L_{n-1}\cdot D\cdot
U_{n-1}U_{n-2}\cdots U_{1}$ (64)
$\displaystyle=\underline{D_{1}L_{1}}L_{2}\cdots L_{n-1}\cdot D\cdot
U_{n-1}U_{n-2}\cdots U_{1},\quad\mbox{since }D_{1}=I,$
$\displaystyle=\overline{L}_{1}\underline{D_{2}L_{2}}\cdots L_{n-1}\cdot
D\cdot U_{n-1}U_{n-2}\cdots U_{1},\quad\mbox{using
}D_{1}L_{1}=\overline{L}_{1}D_{2}\mbox{ as in \eqref{DLLD},}$
$\displaystyle=\ldots$
$\displaystyle=\overline{L}_{1}\overline{L}_{2}\cdots\underline{D_{n-1}L_{n-1}}\cdot
D\cdot U_{n-1}U_{n-2}\cdots U_{1}$
$\displaystyle=\overline{L}_{1}\overline{L}_{2}\cdots\overline{L}_{n-1}\cdot(D_{n}D)\cdot
U_{n-1}U_{n-2}\cdots U_{1}$

(the factors being transformed on each step using (63) are underlined).

The procedure is analogous for the upper bidiagonal factors.

Thus, if (64) is a singularity-free bidiagonal decomposition of the TN matrix
$A$, then

$A=\overline{L}_{1}\overline{L}_{2}\cdots\overline{L}_{n-1}\cdot\overline{D}\cdot\overline{U}_{n-1}\overline{U}_{n-2}\cdots\overline{U}_{1}$

is a singularity-free idiagonal decomposition of $A$ with all $(n,n)$ entries
of $\overline{L}_{i}$ and $\overline{U}_{i}$, $i=1,2,\ldots,n$, equal to 1.

The matrices $\overline{L}_{i},\overline{U}_{i},i=1,2,\ldots,n-1$, and
$\overline{D}$ differ from the matrices $L_{i},U_{i}$ and $D$, respectively,
only in

$(\overline{L}_{i})_{n,n}=(\overline{U}_{i})_{n,n}=1$

and

$\displaystyle(\overline{L}_{i})_{n,n-1}$
$\displaystyle=(L_{i})_{n,n-1}\prod_{k=1}^{i-1}(L_{k})_{n,n},$
$\displaystyle(\overline{U}_{i})_{n,n-1}$
$\displaystyle=(U_{i})_{n,n-1}\prod_{k=1}^{i-1}(U_{k})_{n,n},$
$\displaystyle\overline{D}_{nn}$
$\displaystyle=D_{nn}\prod_{k=1}^{n-1}(L_{k})_{n,n}\prod_{k=1}^{n-1}(U_{k})_{n,n}.$

The routine `STNFixBottomRightOfBD` in our package STNTool [21] implements the
technique from this section.

For example, by implementing the technique from this section to the
singularity-free bidiagonal decomposition of our $3\times 3$ example (40) we
obtain a new singularity-free bidiagonal decomposition with all $(n,n)$
entries of all bidiagonal factors now equal to 1:

$\displaystyle\left[\begin{array}[]{ccc}1&x&x^{2}\\\ 1&y&y^{2}\\\
1&z&z^{2}\end{array}\right]$ $\displaystyle=\left[\begin{array}[]{ccc}1\\\
&1\\\ &1&1\end{array}\right]\left[\begin{array}[]{ccc}1\\\ 1&y-x\\\
&1&1\end{array}\right]\left[\begin{array}[]{ccc}1\\\ &1\\\
&&(z-y)(z-x)\end{array}\right]$ (77) $\displaystyle\phantom{=}\mbox{
}\;\times\left[\begin{array}[]{ccc}1&x\\\ &1&y\\\
&&1\end{array}\right]\left[\begin{array}[]{ccc}1\\\ &1&x\\\
&&1\end{array}\right].$ (84)

##  8 Numerical accuracy

In the standard “$1+\delta$” model of floating point arithmetic [26], to which
the IEEE 754 double precision arithmetic [27] conforms, the result of any
floating point calculation is assumed to satisfy

$\mbox{fl}(a\odot b)=(a\odot b)(1+\delta),$ (85)

where $\odot\in\\{+,-\times,/\\},\,|\delta|\leq\varepsilon$, and $\varepsilon$
is tiny and is called machine precision.

For a computed quantity, $\hat{x}$ to have high relative accuracy, it means
that it satisfies an error bound with its true counterpart, $x$

$|\hat{x}-x|\leq\theta|x|,$

where $\theta$ is a modest multiple of $\varepsilon$. In other words, the sign
and most significant digits of $x$ must be correct. In particular, if $x=0$,
it must be computed exactly.

The above model directly implies that the accuracy in numerical calculations
is lost due to one phenomenon only, known as subtractive cancellation [28]. It
occurs when a subtraction of previously rounded off quantities results in the
loss of significant digits. Multiplication, division, and addition of same-
sign quantities preserve the relative accuracy. Subtraction of initial data
such as the nodes $x_{i}$ is also fine, since initial data is assumed to be
exact: (85) tells us the result of that subtraction is computed to high
relative accuracy. The only subtractions in any singularity-free bidiagonal
decomposition presented in this paper are between exact initial data: either
in the form $x_{i}-x_{j}$ between the nodes $x_{i}$ of the Vandermonde-type
matrix as in (54) or $1-x_{i}$ between the exact double precision floating
point number $1$ and a node $x_{i}$ as in (87), (89), and (90).

Detailed error analyses for the ordinary bidiagonal decompositions of all
nonsingular TN Vandermonde-type matrices presented in this paper have already
been performed in the corresponding papers (e.g., [14] for the Bernstein-
Vandermonde matrices, etc.). All those decompositions are computable to high
relative accuracy componentwise.

The new singularity-free bidiagonal decompositions inherit the same
componentwise error bounds and are thus also computable to high relative
accuracy: the offdiagonal entries in the bidiagonal factors $L_{i}$ and
$U_{i}$ as well as the diagonal entries of the diagonal factor $D$ in the
singularity-free bidiagonal decomposition (44) are the same as the
corresponding entries in the ordinary bidiagonal decomposition (41) except for
the factors $v_{ij}$ from (46). The diagonal entries in the bidiagonal factors
$L_{i}$ and $U_{i}$ in (44) are either equal to 1 or are of the form
$x_{i}-x_{j}$, the latter computable with relative error bounded by the
machine precision, $\varepsilon$, per (85).

##  9 Numerical experiments

We performed extensive numerical tests to verify the correctness of the
formulas we derived in this paper as well as their inherent accuracy in
numerical computations. We present one illustrative example in Figure 1. Using
the formulas in section 5.2, we computed the eigenvalues of the $24\times 24$
$q$-Bernstein-Vandermonde matrix $A$ with parameter $q=0.1$ and nodes

$0.1,0.2,0.2,0.2,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,$
$0.5,0.6,0.7,0.7,0.7,0.7,0.7,0.7,0.8,0.9.$

Figure 1: The eigenvalues of a $24\times 24$ $q$-Bernstein-Vandermonde matrix
as computed by the conventional eigenvalue algorithms of LAPACK as implemented
by MATLAB’s eig and also by forming the accurate singularity-free bidiagonal
decomposition as in section 5 and then using the accurate eigenvalue algorithm
of [8].

A $q$-Bernstein-Vandermonde matrix with distinct nodes is nonsingular [10],
thus the 17 rows corresponding to distinct nodes are linearly independent. The
repeated nodes ($0.2$, three times, and $0.7$, six times) correspond to
repeated rows in the matrix and thus the rank of the matrix is 17.

We started with the singularity-free bidiagonal decomposition of $A$, computed
using the formulas in section 5.2, and implemented in our software package
STNTool [21] as the routine `STNBDqBernsteinVandermonde`. We computed the
eigenvalues using the algorithm `STNEigenValues` [8]. The matrix $A$ is never
explicitly formed or needed in this part of the computation, which is not only
accurate (as we see below), but also very efficient: The singularity-free
bidiagonal decomposition takes $O(n^{2})$ time and the eigenvalues take an
additional $O(n^{3})$ time [8].

For comparison, in double precision floating point arithmetic [27], we formed
$A$ explicitly and computed its eigenvalues using the conventional eigenvalue
algorithm of LAPACK [18] (as implemented by `eig` in MATLAB [29]). As
expected, only the largest eigenvalues of $A$ are computed accurately by
`eig`. Further, since $A$ is unsymmetric and its largest eigenvalue is about
1, the complete loss of relative accuracy occurs in all eigenvalues smaller
than about $10^{-12}$, a value much larger machine precision (about
$10^{-16}$). The zero eigenvalues are lost to roundoff. This behavior is fully
expected and justified for the general, structure-ignoring algorithms of
LAPACK. It also underscores the utility of developing special algorithms for
structured matrices, which deliver results to high relative accuracy with the
same $O(n^{3})$ efficiency.

For verification, we formed $A$, then computed its eigenvalues, in 60 decimal
digit arithmetic using the software package Mathematica. All nonzero
eigenvalues computed by `STNEigenValues` agreed with those computed by
Mathematica to at least $14$ significant decimal digits. No amount of extra
precision will allow us to reliably compute the zero eigenvalues using the
conventional eigenvalue algorithms of LAPACK, thus the only reason we know the
algorithm from [8] computed the correct number of zero eigenvalues ($7$) is
because we know that the rank of the matrix is 17.333Since the graph is log-
scale, the zero eigenvalues are not depicted and the negative and complex
eigenvalues returned by the conventional algorithm are displayed by their
absolute values.

MATLAB implementations of the algorithms described in this paper are available
online [21].

##  10 Appendix

Here we present the explicit formulas for the singularity-free bidiagonal
decompositions of the following Vandermonde-type matrices:

  * 1.

$h$-Bernstein-Vandermonde

  * 2.

Lupaş

  * 3.

rational Bernstein-Vandermonde, and

  * 4.

Cauchy-Vandermonde matrices with one multiple pole.

This is not a complete list of Vandermonde-type matrices as this is an active
area of research, but since we derived these formulas before realizing that
all matrices of the Vandermonde-type follow the same pattern, we share them
here and have implemented them in software [21].

###  10.1 Lupaş matrices

The Lupaş $q$-analogues of the Bernstein functions of degree $n$ for $0<q\leq
1$ are defined as [30]

$l_{i,q}^{n}(x)=\frac{a_{i,q}^{n}(x)}{w_{q}^{n}(x)},\quad x\in[0,1],\quad
i=0,1,\ldots,n,$

where

$\displaystyle a_{i,q}^{n}(x)$
$\displaystyle=\left[\genfrac{}{}{0.0pt}{0}{n}{i}\right]q^{\frac{i(i-1)}{2}}x^{i}(1-x)^{n-i},$
$\displaystyle w_{q}^{n}(x)$
$\displaystyle=\sum_{i=0}^{n}a_{i,q}^{n}(x)=\prod_{i=2}^{n}(1-x+q^{i-1}x).$

For case $q=1$ the Lupaş $q$-analogues of the Bernstein functions coincide
with the Bernstein polynomials.

An $n\times n$ Lupaş matrix $A$ with nodes $x_{1},x_{2},\ldots,x_{n}$ is
defined in [11] as

$A=\left[l^{n-1}_{j-1,q}(x_{i})\right]_{i,j=1}^{n}.$

For $0\leq x_{1}\leq x_{2}\leq\cdots\leq x_{n}<1$, the matrix $A$ is TN [11,
Thm. 2.1]. It is also a Vandermonde-type matrix on $[0,1)$ with a singularity-
free bidiagonal decomposition $[B,C]=\mathcal{SBD}(A)$ defined as in (52)–(55)
with parameters

$\displaystyle s_{ij}$
$\displaystyle=\frac{(1-x_{i})^{n-j}(1-x_{i-j-1})}{(1-x_{i-1})^{n+1-j}}\prod_{k=1}^{n-2}\frac{1-x_{i-1}+q^{k}x_{i-1}}{1-x_{i}+q^{k}x_{i}},$
$\displaystyle m_{ji}$
$\displaystyle=\frac{[n-i+1]q^{i-2}x_{j}}{[i-1](1-x_{j})},$ for $1\leq j<i\leq
n$ and $\displaystyle s_{ii}$
$\displaystyle=\left[\genfrac{}{}{0.0pt}{0}{n-1}{i-1}\right]\frac{q^{\frac{(i-1)(i-2)}{2}}(1-x_{i})^{n-i}}{\prod_{k=2}^{n-1}(1-x_{i}+q^{k-1}x_{i})\prod_{k=0}^{i-1}(1-x_{k})}$

for $1\leq i\leq n$.

The singularity-free bidiagonal decomposition of a Lupaş matrix can be
computed to high relative accuracy in $O(n^{2})$ time using the routine
`STNBDLupas` in our package STNTool [21].

###  10.2 $h$-Bernstein-Vandermonde matrices

These matrices are generalization of the Bernstein-Vandermonde matrices [12,
14], which are the $h=0$ case. The nonsingular case was developed in [15].

The $h$-Bernstein polynomials of degree $n$ for a real parameter $h\geq 0$ are
defined in [15] as

$b_{i,h}^{n}(x)=\left(\genfrac{}{}{0.0pt}{0}{n}{i}\right)\frac{\prod_{k=0}^{i-1}(x+kh)\prod_{k=0}^{n-i-1}(1-x+kh)}{\prod_{k=0}^{n-1}(1+kh)},\;x\in[0,1],\;0\leq
i\leq n.$ (86)

An $n\times n$ $h$-Bernstein-Vandermonde matrix with nodes
$x_{1},x_{2},\ldots,x_{n}$ is defined as

$A=\left[b^{n-1}_{j-1,h}(x_{i})\right]_{i,j=1}^{n}.$

For $0\leq x_{1}\leq x_{2}\leq\cdots\leq x_{n}<1$, $A$ is a TN Vandermonde-
type matrix [15] on $[0,1)$ and has a singularity-free bidiagonal
decomposition $[B,C]=\mathcal{SBD}(A)$ defined as in (52)–(55) with parameters

$\displaystyle s_{ij}$
$\displaystyle=\frac{1-x_{i-j}+(n-j)h}{1-x_{i-1}+(n-j)h}\prod_{k=0}^{n-j-1}\frac{1-x_{i}+kh}{1-x_{i-1}+kh},$
(87) $\displaystyle m_{ji}$
$\displaystyle=\frac{n-i+1}{i-1}\cdot\frac{(x_{j}+(i-j-1)h)\prod_{k=1}^{j-1}(1-x_{k}+(n-i+1)h)}{\prod_{k=1}^{j}(1-x_{k}+(n-i)h)},$
for $1\leq j<i\leq n$ and $\displaystyle s_{ii}$
$\displaystyle=\left(\genfrac{}{}{0.0pt}{0}{n-1}{i-1}\right)\frac{\prod_{k=0}^{n-i-1}(1-x_{i}+kh)}{\prod_{k=1}^{n-i-1}(1+kh)\prod_{k=1}^{i-1}(1-x_{k}+(n-i)h)}$

for $1\leq i\leq n$.

The singularity-free bidiagonal decomposition of an $h$-Bernstein-Vandermonde
matrix can be computed to high relative accuracy in $O(n^{2})$ time using the
routine `STNBDhBernsteinVandermonde` in our package STNTool [21].

###  10.3 Rational Bernstein-Vandermonde matrices

A rational Bernstein-Vandermonde matrix with nodes $x_{1},x_{2},\ldots,x_{n}$
and positive weights $w_{1},w_{2},\ldots,w_{n}$ is defined as [9]:

$R=\left[r_{j-1}^{n-1}(x_{i})\right]_{i,j=1}^{n},$

where

$r_{j-1}^{n-1}(x)=\frac{w_{j}\,b_{j-1}^{n-1}(x)}{W(x)},\quad x\in[0,1],\quad
j\in\\{1,2,\ldots,n\\},$ (88)

and $W(x)=\sum_{j=1}^{n}w_{j}\,b_{j-1}^{n-1}(x)$. The functions $b_{j}^{n}(x)$
are the Bernstein polynomials of degree $n$. They equal the $q$-Bernstein
polynomials (57) for $q=1$ and the $h$-Bernstein polynomials (86) for $h=0$,
i.e.,

$b_{i}^{n}(x)=\left(\genfrac{}{}{0.0pt}{0}{n}{i}\right)x^{i}(1-x)^{n-i},\quad
i\in\\{0,1,\ldots,n\\},\quad x\in[0,1].$

For $0<x_{1}\leq x_{2}\leq\cdots\leq x_{n}<1$, the matrix $R$ is TN [9, Thm.
3.1], which is of Vandermonde-type on $(0,1)$ and has a singularity-free
bidiagonal decomposition $[B,C]=\mathcal{SBD}(R)$ defined as in (52)–(55) with
parameters

$\displaystyle s_{ij}$
$\displaystyle=\frac{W(x_{i-1})(1-x_{i})^{n-j}(1-x_{i-j})}{W(x_{i})(1-x_{i-1})^{n-j+1}},$
(89) $\displaystyle m_{ji}$
$\displaystyle=\frac{w_{i}}{w_{i-1}}\cdot\frac{n-i+1}{i-1}\cdot\frac{x_{j}}{1-x_{j}},$
(90) for $1\leq j<i\leq n$ and $\displaystyle s_{ii}$
$\displaystyle=\left(\genfrac{}{}{0.0pt}{0}{n-1}{i-1}\right)\frac{w_{i}(1-x_{i})^{n-i}}{W(x_{i})\prod_{k=1}^{i-1}(1-x_{k})}$

for $1\leq i\leq n$.

While we like the above expressions for their symmetry, an alternative
singularity-free bidiagonal decomposition of $R$ can be obtained by
recognizing, from (58) and (88), that $R=W_{1}^{-1}AW_{2},$ where $A$ is the
$q$-Bernstein-Vandermonde matrix for $q=1$, and

$\displaystyle W_{1}$ $\displaystyle={\rm
diag}(W(x_{1}),W(x_{2}),\ldots,W(x_{n}))$ $\displaystyle W_{2}$
$\displaystyle={\rm diag}(w_{1},w_{2},\ldots,w_{n})$

are positive diagonal matrices. Thus if $\mathcal{SBD}(A)$ is given by (44),
then $R$ has a singularity-free bidiagonal decomposition

$R=(W_{1}^{-1}L_{1})L_{2}\cdots L_{n-1}DU_{n-1}U_{n-2}\cdots(U_{1}W_{2}),$
(91)

where the first and last factors in the parentheses, $W_{1}^{-1}L_{1}$ and
$U_{1}W_{2}$,, are nonnegative bidiagonal matrices with the same nonzero
pattern as $L_{1}$ and $U_{1}$, respectively. Thus (91) is another bidiagonal
decomposition of $R$ and an example of how the nonuniqueness of the
singularity-free bidiagonal decomposition (44) can play out.

The singularity-free bidiagonal decomposition of a rational Bernstein-
Vandermonde matrix can be computed to high relative accuracy in $O(n^{2})$
time using the routine `STNBDRationalBernsteinVandermonde` in our package
STNTool [21].

###  10.4 Cauchy–Vandermonde matrices with one multiple pole

A Cauchy-Vandermonde matrix with nodes $x_{1},\ldots,x_{n}$ and one pole $d$
of multiplicity $s\geq 1$ is defined as [31]

$A=\left[\begin{array}[]{ccccccc}\frac{1}{(x_{1}+d)^{s}}&\cdots&\frac{1}{x_{1}+d}&1&x_{1}&\cdots&x_{1}^{n-s-1}\\\
\frac{1}{(x_{2}+d)^{s}}&\cdots&\frac{1}{x_{2}+d}&1&x_{2}&\cdots&x_{2}^{n-s-1}\\\
\vdots&&\vdots&\vdots&\vdots&&\vdots\\\
\frac{1}{(x_{n}+d)^{s}}&\cdots&\frac{1}{x_{n}+d}&1&x_{n}&\cdots&x_{n}^{n-s-1}\\\
\end{array}\right]$

and is totally nonnegative when $0\leq x_{1}\leq x_{2}\leq\cdots\leq x_{n}$
and $d>0$. It is a Vandermonde-type matrix on $[0,\infty)$ for $d>0$ and has a
decomposition $[B,C]=\mathcal{SBD}(A)$ with parameters

$\displaystyle s_{ij}$
$\displaystyle=\left(\frac{x_{i-1}+d}{x_{i}+d}\right)^{s},$ $\displaystyle
m_{ji}$ $\displaystyle=\left\\{\begin{array}[]{ll}x_{j}+d,&\hbox{$0<i-j\leq
s$}\\\ x_{j},&\hbox{$i-j>s$}\end{array}\right.,$ for $1\leq j<i\leq n$ and
$\displaystyle s_{ii}$ $\displaystyle=\frac{1}{(x_{i}+d)^{s}}$

for $1\leq i\leq n$.

The singularity-free bidiagonal decomposition of a Cauchy-Vandermonde matrix
with one multiple pole can be computed to high relative accuracy in $O(n^{2})$
time using the routine `STNBDCauchyVandermonde1pole` in our package STNTool
[21].

###  10.5 Other Vandermonde-type matrices

The method presented in this section can be used analogously to derive
singularity-free bidiagonal decompositions of other Vandermonde-type matrices,
for example:

  * 1.

Said-Ball Vandermonde [13, 16],

  * 2.

Rational Said-Ball Vandermonde [9],

  * 3.

Collocation and Wronskian matrices of Jacobi polynomials [32].

This is likely an incomplete list as this currently appears to be an area of
active research.

##  11 Acknowledgements

We thank the anonymous referees for the very careful reading of our manuscript
as well as their very constructive comments and suggestions, which greatly
improved the presentation and clarity of the paper and, in particular, for
pointing out the significance of the result now formulated in Corollary 4.1.
We also thank Slobodan Simić for his help with the proof of that corollary.

P. Koev thanks the University of California – Berkeley for the kind
hospitality during his sabbatical leave in 2019.

## References

  * [1] T. Ando, Totally positive matrices, Linear Algebra Appl. 90 (1987) 165–219.
  * [2] S. M. Fallat, C. R. Johnson, Totally nonnegative matrices, Princeton Series in Applied Mathematics, Princeton University Press, Princeton, NJ, 2011.
  * [3] S. Karlin, Total Positivity. Vol. I, Stanford University Press, Stanford, CA, 1968.
  * [4] S. M. Fallat, Bidiagonal factorizations of totally nonnegative matrices, Amer. Math. Monthly 108 (8) (2001) 697–712.
  * [5] M. Gasca, J. M. Peña, Total positivity and Neville elimination, Linear Algebra Appl. 165 (1992) 25–44.
  * [6] P. Koev, Accurate eigenvalues and SVDs of totally nonnegative matrices, SIAM J. Matrix Anal. Appl. 27 (1) (2005) 1–23.
  * [7] P. Koev, Accurate computations with totally nonnegative matrices, SIAM J. Matrix Anal. Appl. 29 (2007) 731–751.
  * [8] P. Koev, Accurate eigenvalues and zero Jordan blocks of (singular) totally nonnegative matrices, Numer. Math. 141 (2019) 693–713.
  * [9] J. Delgado, J. M. Peña, Accurate computations with collocation matrices of rational bases, Appl. Math. Comput. 219 (9) (2013) 4354–4364.
  * [10] J. Delgado, J. M. Peña, Accurate computations with collocation matrices of q-Bernstein polynomials, SIAM J. Matrix Anal. Appl. 36 (2) (2015) 880–893.
  * [11] J. Delgado, J. M. Peña, Accurate computations with Lupaş matrices, Appl. Math. Comput. 303 (2017) 171–177.
  * [12] A. Marco, J. J. Martínez, A fast and accurate algorithm for solving Bernstein-Vandermonde linear systems, Linear Algebra Appl. 422 (2-3) (2007) 616–628.
  * [13] A. Marco, J. J. Martínez, Accurate computations with Said-Ball-Vandermonde matrices, Linear Algebra Appl. 432 (11) (2010) 2894–2908.
  * [14] A. Marco, J. J. Martínez, Accurate computations with totally positive Bernstein-Vandermonde matrices, Electron. J. Linear Algebra 26 (2013) 357–380.
  * [15] A. Marco, J. J. Martínez, R. Viaña, Accurate bidiagonal decomposition of totally positive h-Bernstein-Vandermonde matrices and applications, Linear Algebra Appl. 579 (2019) 320–335.
  * [16] A. Marco, J. J. Martínez, Bidiagonal decomposition of rectangular totally positive Said-Ball-Vandermonde matrices: error analysis, perturbation theory and applications, Linear Algebra Appl. 495 (2016) 90–107.
  * [17] L. Han, Y. Chu, Z. Qiu, Generalized Bézier curves and surfaces based on Lupaş q-analogue of Bernstein operator, Journal of Computational and Applied Mathematics 261 (2014) 352–363.
  * [18] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users’ Guide, Third Edition, Software Environ. Tools 9, SIAM, Philadelphia, 1999.
  * [19] J. W. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia, 1997.
  * [20] M. J. Ablowitz, A. S. Fokas, Complex variables: introduction and applications, 2nd Edition, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, 2003.
  * [21] P. Koev, `http://www.math.sjsu.edu/~koev`.
  * [22] G. M. Phillips, Bernstein polynomials based on the $q$-integers, Ann. Numer. Math. 4 (1-4) (1997) 511–518.
  * [23] J. Delgado, H. Orera, J. M. Peña, Accurate computations with Laguerre matrices, Numer. Linear Algebra Appl. 26 (1) (2019) e2217, 10.
  * [24] J. Delgado, H. Orera, J. M. Peña, Accurate algorithms for Bessel matrices, J. Sci. Comput. 80 (2) (2019) 1264–1278.
  * [25] E. Mainar, J. M. Peña, B. Rubio, Accurate computations with Wronskian matrices, Calcolo 58 (1) (2021) Paper No. 1, 15.
  * [26] N. J. Higham, Accuracy and Stability of Numerical Algorithms, Second Edition, SIAM, Philadelphia, 2002.
  * [27] ANSI/IEEE, New York, IEEE Standard for Binary Floating Point Arithmetic, Std 754-1985 Edition (1985).
  * [28] J. Demmel, Accurate singular value decompositions of structured matrices, SIAM J. Matrix Anal. Appl. 21 (2) (1999) 562–580.
  * [29] The MathWorks, Inc., Natick, MA, MATLAB Reference Guide (1992).
  * [30] A. Lupaş, A $q$-analogue of the Bernstein operator, in: Seminar on Numerical and Statistical Calculus (Cluj-Napoca, 1987), Vol. 87 of Preprint, Univ. “Babeş-Bolyai”, Cluj-Napoca, 1987, pp. 85–92.
  * [31] J. Martínez, J. M. Peña, Factorizations of Cauchy-Vandermonde matrices with one multiple pole, in: Recent research on pure and applied algebra, Nova Science Publishers, Hauppauge, NY, 2003, pp. 85–95.
  * [32] E. Mainar, J. M. Peña, B. Rubio, Accurate computations with collocation and Wronskian matrices of Jacobi polynomials, J. Sci. Comput. 87 (3) (2021) Paper No. 77, 30.


# On problems of the Divergence of logarithmic means of Fourier series

Ushangi Goginava  U. Goginava, Department of Mathematical Sciences
United Arab Emirates University, P.O. Box No. 15551
Al Ain, Abu Dhabi, UAE<EMAIL_ADDRESS><EMAIL_ADDRESS>and
Farrukh Mukhamedov  F. Mukhamedov, New Uzbekistan University, 54, Mustaqillik
ave.
100007, Tashkent, Uzbekistan; Central Asian University, Tashkent 111221,
Uzbekistan; United Arab Emirates University, P.O. Box No. 15551
Al Ain, Abu Dhabi, UAE<EMAIL_ADDRESS>
###### Abstract.

The main objective of the present paper is to solve the several long standing
divergence problems related to the logarithmic means of Fourier series in
context of general orthonormal systems.

00footnotetext: 2020 Mathematics Subject Classification: 42C10, 40A05, 40A30,
42A24, 42A55. Key words and phrases: Almost Everywhere Divergence, Logarithmic
Means, Fourier Series.

##  1\. Introduction

Let $\psi:[0,\infty)\rightarrow[0,\infty)$ be a non-decreasing function. The
class of all measurable functions $f:\mathbb{I}\rightarrow\mathbb{R}$,
$\mathbb{I}:=[0,1)$ for which $\int_{\mathbb{I}}\psi(|f|)<\infty$ will be
denoted by $\psi\left(L\right)\left(\mathbb{I}\right)$. If
$\psi\left(u\right)=u$ then $\psi\left(L\right)\left(\mathbb{I}\right)$
coincides with the standard $L_{1}$-space, i.e.
$\psi\left(L\right)\left(\mathbb{I}\right):=L_{1}\left(\mathbb{I}\right)$.

In the sequel, by $\Psi$ we denote the set of all non-decreasing functions
$\psi:[0,\infty)\rightarrow[0,\infty)$ such that
$\liminf\limits_{u\rightarrow\infty}\psi(u)/u>0$. Note that for every set
$E\subset\mathbb{R}$ with finite measure and every function $\psi\in\Psi$, the
inclusion $\psi(L)(E)\subset L_{1}(E)$ is valid.

It is well-known that Kolmogoroff in [22] solved the Lusin’s problem by
constructing an integrable function $f\in L_{1}(\mathbb{I})$ whose Fourier
trigonometric series diverges almost everywhere. Moreover, in his work [23],
Kolmogorov provided an example of a function $f\in L_{1}(\mathbb{I})$ where
the divergence occurs everywhere. Similar outcomes were demonstrated by Stein
for the Walsh-Paley system, as detailed in [32]. These results can be
summarized as follows.

######  Theorem KS.

Assume that $S_{n}(f,x)$ is the partial sum of Fourier series of a function
$f$ with respect to either the trigonometric or the Walsh-Paley system. Then
there is $f\in L_{1}\left(\mathbb{I}\right)$ such that

(1) $\sup_{n}|S_{n}(f,x)|=\infty,\ \ \ \text{for all}\ x\in\mathbb{I}.$

Let $\Phi={\varphi_{n}(x)}$ denote an orthonormal system, satisfying
$|\varphi_{n}(x)|\leq M$ for all $x\in\mathbb{I}$ and $n\in\mathbb{N}$.
Bochkarev [2] has established the existence of $f\in
L_{1}\left(\mathbb{I}\right)$ such that the inequality (1) holds for the
partial sums of the Fourier series of $f$ concerning the system $\Phi$ on a
certain set $E\subset\mathbb{I}$ with $|E|>0$. It is important to note that
Kazarian [21] demonstrated that the result by Bochkarev does not hold for the
set $E$ with $|E|=1$.

The exploration into the divergence of Fourier series remains an ongoing
pursuit, with the challenge of identifying the precise class of functions
exhibiting divergence. In this context, we highlight certain established
findings in this direction.

In [25] Konyagin established, for the trigonometric system, the existence of a
function

$f\in\psi(L)\left(\mathbb{I}\right)$

such that $\sup_{k\in\mathbb{N}}\left|S_{k}(f)(x)\right|=\infty$ holds for
every $x\in\mathbb{I}$, where $\psi(u)=o(u\sqrt{\log u/\log\log u})$.

If $\psi(u)=o(u\sqrt{\log u})(u\rightarrow\infty)$, then Bochkarev [3]
investigated the Walsh-Paley system and established the existence of a
function $f\in\psi(L)\left(\mathbb{I}\right)$ such that

$\sup_{k\in\mathbb{N}}\left|S_{k}(f)(x)\right|=\infty$

holds for every $x\in\mathbb{I}$.

We recall that the logarithmic means of a function is defined by

(2) $L_{n}(f,x)=\frac{1}{\ell_{n}}\sum_{j=0}^{n-1}\frac{S_{j}(f,x)}{n-j},$

where $\ell_{n}:=\sum\limits_{k=1}^{n}\left(1/k\right)$. It is well-known that
the convergence of the logarithmic means is weaker than the convergence of the
partial sums. Consequently, such types of means have been applied in the
Fourier series. In 1992, the convergence of logarithmic means of Fourier
series was posed as a problem in [30]. While this problem has been
successfully addressed for the Walsh system [11, 12], the question of almost
everywhere convergence/divergence for the trigonometric or general orthonormal
system remains open.

The principal objective of this paper is to resolve the aforementioned issues
concerning both trigonometric and general orthonormal systems. To achieve
this, we have derived a concise and simplified result (refer to Theorem 6),
enabling us to establish our key findings. Let us formulate our first result
of this paper.

######  Theorem 1.

The next statements hold true:

  * (i)

Let $\psi(u)=o(u\sqrt{\log u/\log\log u})$ then there is
$f\in\psi(L)\left(\mathbb{I}\right)$ such that

(3) $\sup_{n}|L_{n}(f,x)|=\infty,\ \ \ \text{for all}\ x\in\mathbb{I},$

for the trigonometric system;

  * (ii)

Let $\psi(u)=o(u\sqrt{\log u})$ then there is
$f\in\psi(L)\left(\mathbb{I}\right)$ such that (3) holds for the Walsh-Paley
system;

  * (iii)

Let $\Phi=\\{\varphi_{n}(x)\\}$ be an orthonormal system such that
$|\varphi_{n}(x)|\leq M$ for all $x\in\mathbb{I}$, $n\in\mathbb{N}$. Then
there is $f\in L_{1}\left(\mathbb{I}\right)$ and a set $E\subset\mathbb{I}$
with $|E|>0$ such that

$\sup_{n}|L_{n}(f,x)|=\infty,\ \ \ \text{for all}\ x\in E.$

It is important to emphasize that part (ii) of Theorem 1 had been previously
demonstrated in [12], where the divergence of logarithmic means was
established only on a set of positive measure. However, Theorem 1 not only
enhances the previously proven result in [12] but also offers a more concise
proof for it.

On a different note, Gosselin [15] established that for every increasing
sequence of natural numbers $\left(n_{k}\right)$, there exists a function
$f\in L_{1}(\mathbb{I})$ such that

(4) $\sup_{k\in\mathbb{N}}\left|S_{n_{k}}(f,x)\right|=\infty$

for almost every $x\in\mathbb{I}$. Conversely, Totik [34] found a function
$f\in L_{1}(\mathbb{I})$ for which (4) holds for every $x\in\mathbb{I}$,
considering any increasing sequence of natural numbers $\left(n_{k}\right)$.

It turns out that establishing an analogue of the results by Gosselin and
Totik remains an open problem, intricately linked with Konyagin’s problem
[24]. While a partial solution has been provided in [14]. Let us recall some
necessary notions to formulate it.

The _binary coefficients_ of a number $n\in\mathbb{N}$ will be denoted by
$\varepsilon_{j}(n)$ $(j\in\mathbb{N}_{0})$, i.e.,
$\varepsilon_{j}(n)\in\\{0,1\\}$ for every $j\in\mathbb{N}_{0}$ and
$n=\sum_{j=0}^{\infty}\varepsilon_{j}(n)2^{j}$. Note that
$\varepsilon_{j}(n)=0$ if $j$ is sufficiently large. For $n\in\mathbb{N}$, let
us denote

$V(n)=\varepsilon_{0}(n)+\sum_{j=1}^{\infty}|\varepsilon_{j}(n)-\varepsilon_{j-1}(n)|.$

The quantity $V(n)$ is called the _variation_ of the number $n$.

Let us define the _spectrum_ of a number $n\in\mathbb{N}$ as follows

$\text{Sp}(n)=\\{j\in\mathbb{N}:\varepsilon_{j}(n)=1\\}.$

Let $(n_{k})$ be an increasing sequences of natural numbers. We will call
$(n_{k})$ a sequence with _nested spectrums_ if

$\text{Sp}(n_{k+1})\cap[0,\max\text{Sp}(n_{k})]=\text{Sp}(n_{k})\;\;\;\text{for
every}\;\;\;k\in\mathbb{N}.$

######  Theorem GO.

[14] Let $(n_{k})$ be a sequence of natural numbers with nested spectrums
having unbounded variation. Then there exists a function $f\in
L_{1}\left(\mathbb{I}\right)$ such that

$\sup_{k\in\mathbb{N}}|S_{n_{k}}(f,x)|=\infty$

for every $x\in\mathbb{I}.$

Next our result improves the above mentioned Theorems for the logarithmic
means.

######  Theorem 2.

The next statements hold true:

  * (i)

for every increasing sequence of natural numbers $\left(n_{k}\right)$ there
exists a function $f\in L_{1}(\mathbb{I})$ such that

(5) $\sup_{N}\bigg{|}\frac{1}{\ln
N}\sum_{k=0}^{N-1}\frac{S_{n_{k}}(f,x)}{N-k}\bigg{|}=\infty,\ \ \ \text{for
all}\ x\in\mathbb{I}.$

for the trigonometric system;

  * (ii)

Let $(n_{k})$ be nested sequence which has unbounded variation. Then there is
$f\in L_{1}(\mathbb{I})$ such that (5) holds for the Walsh-Paley system.

######  Remark 1.

Recently, Gát [10] has provided necessary and sufficient conditions on
subsequence $(n_{k})$, for the trigonometric system, such that

(6) $\frac{1}{N}\sum_{k=0}^{N-1}S_{n_{k}}(f,x)\rightarrow f\ \ \text{a.e.}$

for every $f\in L_{1}(\mathbb{I})$. However, from Theorem 1 we infer that the
arithmetic means in (6) cannot be replaced by the logarithmic ones.

Now, let us consider any lacunary sequence $(n_{k})$. From Theorem 2 one can
ask: find an optimal class of functions for which the sequence

$\frac{1}{\ln N}\sum_{k=0}^{N-1}\frac{S_{n_{k}}(f,x)}{N-k}$

converges a.e. to $f$. Our result allows us to provide an affirmative answer
to this question. Namely, we can prove the next result.

######  Theorem 3.

Let $(n_{k})$ be an arbitrary lacunary sequence. Then the space
$L\log^{+}\log^{+}L\log^{+}\log^{+}\log^{+}\log^{+}L\left(\mathbb{I}\right)$
is the optimal Orlicz space for a.e. convergence of

$\frac{1}{l_{N}}\sum_{k=0}^{N-1}\frac{S_{n_{k}}(f,x)}{N-k}\rightarrow f\text{
\ as }n\rightarrow\infty.$

The proof follows immediately from Theorem 6 and the results of Di Plinio [6]
and Lie [29].

On the other hand, Buzdalin [5] provided an affirmative response to Ul’yanov’s
question by demonstrating that for any set $\mathrm{E}\subset\mathbb{I}$ with
$|\mathrm{E}|=0$, there exists a real continuous function
$f\in\mathrm{C}(\mathbb{I})$ such that its Fourier series with respect to the
trigonometric system diverges unboundedly on $E$. Now, leveraging Theorem 6,
we can derive our subsequent result.

######  Theorem 4.

For any set $\mathrm{E}\subset\mathbb{I}$ with $|\mathrm{E}|=0$ there is a
real function $f\in\mathrm{C}(\mathbb{I})$ such that its logarithmic means (2)
of Fourier series with respect to trigonometric system diverges unboundedly on
$E$.

Let us point out that an analogue of Buzdalin’s theorem remains an open
problem for the Walsh systems. Conversely, for bounded functions, this problem
has been resolved by Bughadze [4], and the divergence of arithmetic means on
sets of measure zero has been explored in [13]. Karagulyan [19] investigated a
more general sequence of operators and studied its divergence for bounded
functions on a measure zero set.

It is noteworthy that the almost everywhere convergence and divergence of the
partial sums of Fourier series have been investigated by many authors,
including Setchkin [31], Katznelson [20], Kahane and Katznelson [17], Konyagin
[28, 27, 26], Do and Lacey [7], and Gát [8, 9].

##  2\. Proof of Main Result

Let $\mathbf{S}:=\left\\{s_{n}:n\in\mathbb{N}\right\\}$ be a sequence of real
numbers. The logarithmic means of the sequence $\mathbf{S}$ are defined as
follows

(7)
$L_{n}\left(\mathbf{S}\right):=\frac{1}{\ell_{n}}\sum\limits_{k=0}^{n-1}\frac{s_{k}}{n-k}.$

where, as before, $\ell_{n}:=\sum\limits_{k=1}^{n}\left(1/k\right)$.

Now, let us consider more general means. Assume that
$\mathbf{T}:=\left(t_{k,n}\right)$ is an infinite matrix. Define the $n$th
matrix transform of the sequence $\mathbf{S}$ by

(8)
$T_{n}\left(\mathbf{S}\right):=\sum_{k=0}^{n}t_{k,n}s_{k}\quad(n\in{\mathbb{N}}).$

We say that $T$ is regular if $T_{n}\left(\mathbf{S}\right)\rightarrow
s\left(n\rightarrow\infty\right)$ whenever $s_{n}\rightarrow
s\left(n\rightarrow\infty\right).$ In order that $T$ be regular it is
neccessary and sufficient that (see [35, Ch. 3] or [16, P. 43-57 ]):

  * (a)

$\sum_{k=0}^{n}t_{k,n}=1;$

  * (b)

$t_{k,n}\rightarrow 0$ for each $k$, when $n\rightarrow\infty;$

  * (c)

$\sup\limits_{n}\sum\limits_{k=0}^{n}\left|t_{k,n}\right|\leq C<\infty$.

We notice that if $t_{n,k}=q_{n-k}/Q_{n}$, where $(q_{n})$ is some sequence,
and $Q_{n}:=\sum\limits_{k=0}^{n}q_{k}$, then $T_{n}\left(\mathbf{S}\right)$
coincidences with the Nörlund means, i.e.

(9)
$T_{n}\left(\mathbf{S}\right)=:N_{n}^{\left(q\right)}\left(\mathbf{S}\right)=\frac{1}{Q_{n}}\sum\limits_{k=0}^{n}q_{n-k}s_{k},$

If $q_{n}>0$ for all $n\geq 0$, then the regularity of the Nörlund means is
equivalent to

(10) $\lim_{n\to\infty}\frac{q_{n}}{Q_{n}}=0.$

The $\left(C,\alpha_{n}\right)$ means (see [18] for more details) of the
sequence $\mathbf{S}:=\left\\{s_{n}:n\in\mathbb{P}\right\\}$ is given by

$\sigma_{n}^{\alpha_{n}}(\mathbf{S})=\frac{1}{A_{n}^{\alpha_{n}}}\sum\limits_{j=0}^{n}A_{n-j}^{\alpha_{n}-1}s_{j},$

where

$A_{n}^{\alpha_{n}}:=\frac{(1+\alpha_{n})\dots(n+\alpha_{n})}{n!}.$

The connections between method $\left(C,\alpha_{n}\right)$ with varying
parameters and Cesaro’s method $\left(C,\alpha\right)$ with constant
parameters were studied by Kaplan [18]. Regarding, the connection between two
different methods $\left(C,\alpha_{n}\right)$ and $\left(C,\beta_{n}\right)$
was established recently by Akhobadze and Zviadadze [1].

In the sequel, we are going to employ next well-known result (see [16, pp.
68-69]).

######  Theorem H.

Let $(q_{n})$ be a sequence such that

(11)
$q_{0}=1,q_{n}>0,\frac{q_{n+1}}{q_{n}}\geq\frac{q_{n}}{q_{n-1}}\left(n>0\right),$

Assume that (10) holds. Denote

$q(x)=\sum_{n=0}^{\infty}q_{n}x^{n}.$

Then

(12) $\frac{1}{q\left(x\right)}=\sum\limits_{n=0}^{\infty}\gamma_{n}x^{n},\ \
\ |x|<1,$

where $\gamma_{0}=1$, $\gamma_{n}<0\left(n>0\right),$ and

(13) $\sum\limits_{n=1}^{\infty}\gamma_{n}\geq-1.$

Using the previous result we are going to prove next auxiliary fact.

######  Theorem 5.

Let $\mathbf{S}:=\left\\{s_{n}:n\in\mathbb{P}\right\\}$ be a sequence of real
numbers such that

$L_{n}\left(\mathbf{S}\right)\rightarrow s\text{ \ as \ }n\rightarrow\infty.$

Then for $\alpha_{n}\in\left(0,1\right)$

$\sigma_{n}^{\alpha_{n}}(\mathbf{S})\rightarrow s\text{ \ \ as \
}n\rightarrow\infty.$

.

###### Proof.

Let $m$ be any fixed number. We first notice that the condition (10) (see
[16]) implies the convergence of the series
$\sum\limits_{n=0}^{\infty}s_{n}x^{n}$ when $\left|x\right|<1$. Moreover, from
$\lim\limits_{n\rightarrow\infty}\frac{A_{n}^{\alpha_{m}-1}}{A_{n-1}^{\alpha_{m}-1}}=1$
we infer the convergence of the series
$\sum\limits_{n=0}^{\infty}A_{n}^{\alpha_{m}-1}x^{n}$ when $\left|x\right|<1$.
Hence, by the Mertens Theorem one concludes that the next power series is
convergent (when $|x|<1$)

(14)
$\displaystyle\left(\sum\limits_{n=0}^{\infty}A_{n}^{\alpha_{m}-1}x^{n}\right)\left(\sum\limits_{n=0}^{\infty}s_{n}x^{n}\right)$
$\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}\left(\sum\limits_{k=0}^{n}A_{n-k}^{\alpha_{m}-1}s_{k}\right)x^{n}$
$\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}A_{n}^{\alpha_{m}}\sigma_{n}^{\alpha_{m}}(\mathbf{S})x^{n}.$

Let the sequence $\left(q_{n}\right)_{n=1}^{\infty}$ satisfy the conditions of
Theorem H. Then by Theorem H together with (14) we obtain

$\displaystyle\sum\limits_{n=0}^{\infty}A_{n}^{\alpha_{m}}\sigma_{n}^{\alpha_{m}}(\mathbf{S})x^{n}$
$\displaystyle=$
$\displaystyle\frac{\left(\sum\limits_{n=0}^{\infty}A_{n}^{\alpha_{m}-1}x^{n}\right)}{\left(\sum\limits_{n=0}^{\infty}q_{n}x^{n}\right)}\left(\sum\limits_{n=0}^{\infty}q_{n}x^{n}\right)\left(\sum\limits_{n=0}^{\infty}s_{n}x^{n}\right)$
$\displaystyle=$
$\displaystyle\left(\sum\limits_{n=0}^{\infty}A_{n}^{\alpha_{m}-1}x^{n}\right)\left(\sum\limits_{n=0}^{\infty}\gamma_{n}x^{n}\right)\sum\limits_{n=0}^{\infty}\left(\sum\limits_{k=0}^{n}q_{n-k}s_{k}\right)x^{n}$
$\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}\left(\sum\limits_{k=0}^{n}A_{n-k}^{\alpha_{m}-1}\gamma_{k}\right)x^{n}\sum\limits_{n=0}^{\infty}Q_{n}N_{n}^{\left(q\right)}\left(\mathbf{S}\right)x^{n}.$

Denote

$b_{n,m}:=\sum\limits_{k=0}^{n}A_{n-k}^{\alpha_{m}-1}\gamma_{k}.$

From (14) it follows that

$\displaystyle\sum\limits_{n=0}^{\infty}A_{n}^{\alpha_{m}}\sigma_{n}^{\alpha_{m}}(\mathbf{S})x^{n}$
$\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}b_{n,m}x^{n}\sum\limits_{n=0}^{\infty}Q_{n}N_{n}^{\left(q\right)}\left(\mathbf{S}\right)x^{n}$
$\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}\left(\sum\limits_{k=0}^{n}b_{n-k,m}Q_{k}N_{k}^{\left(q\right)}\left(\mathbf{S}\right)\right)x^{n}.$

Now, we may assume that $n=m$. Then by (2) one finds

$\sigma_{n}^{\alpha_{n}}(\mathbf{S})=\sum\limits_{k=0}^{n}\frac{b_{n-k,n}Q_{k}}{A_{n}^{\alpha_{n}}}N_{k}^{\left(q\right)}\left(\mathbf{S}\right).$

Therefore, the sequence
$\left(\sigma_{n}^{\alpha_{n}}(\mathbf{S})\right)_{n=1}^{\infty}$ can be
represented as a matrix transformation
$\left(t_{k,n}:=\frac{b_{n-k,n}Q_{k}}{A_{n}^{\alpha_{n}}}\right)$ of the
sequence
$\left(N_{k}^{\left(q\right)}\left(\mathbf{S}\right)\right)_{n=1}^{\infty}$.

Let us establish the regularity of this transformation. It is enough to prove
that conditions $\left(a\right)-\left(c\right)$ hold for the sequence
$t_{k,n}$. Indeed, assume that $\mathbf{S}$ is an identical sequence
$s_{n}=1,n\in\mathbb{N}$. Since
$\sigma_{n}^{\alpha_{n}}(\mathbf{S})=N_{k}^{\left(q\right)}\left(\mathbf{S}\right)=1$
we get $\sum\limits_{k=0}^{n}t_{k,n}=1.$ Hence, condition $(a)$ is satisfied.
Now, we prove the condition $\left(b\right)$. For fixed $k$ we can write

$t_{k,n}=\frac{A_{n-k}^{\alpha_{n}}Q_{k}}{A_{n}^{\alpha_{n}}}\frac{1}{A_{n-k}^{\alpha_{n}}}\sum\limits_{l=0}^{n-k}A_{n-k-l}^{\alpha_{n}-1}\gamma_{l}.$

Due to $\frac{A_{n-k}^{\alpha_{n}}Q_{k}}{A_{n}^{\alpha_{n}}}\leq C_{k}$, and
the regularity of the method $\left(C,\alpha_{n}\right)$, for
$\gamma_{k}\rightarrow 0\left(k\rightarrow\infty\right)$ (see (13)) we infer
that

$\frac{1}{A_{n-k}^{\alpha_{n}}}\sum\limits_{l=0}^{n-k}A_{n-k-l}^{\alpha_{n}-1}\gamma_{l}\rightarrow
0$

as $n\rightarrow\infty$. Hence, the condition $\left(b\right)$ holds. It
remains to prove the condition $(c)$.

Assume that the sequence $\left(q_{n}\right)_{n=1}^{\infty}$ satisfies the
following additional condition

(16)
$\frac{q_{n}}{q_{n-1}}\leq\frac{A_{n}^{\alpha_{n}-1}}{A_{n-1}^{\alpha_{n}-1}}=1-\frac{1-\alpha_{n}}{n}\left(n\in\mathbb{N}\right).$

From (12) it follows that

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\left(\sum\limits_{n=0}^{\infty}\gamma_{n}x^{n}\right)\left(\sum\limits_{n=0}^{\infty}q_{n}x^{n}\right)$
$\displaystyle=$
$\displaystyle\sum\limits_{n=0}^{\infty}\left(\sum\limits_{k=0}^{n}q_{n-k}\gamma_{k}\right)x^{n}.$

Consequently,

$\sum\limits_{k=0}^{n}q_{n-k}\gamma_{k}=0\left(n>0\right).$

Then

$\displaystyle\frac{A_{n}^{\alpha_{n}}}{Q_{k}}t_{k,n}$ $\displaystyle=$
$\displaystyle\sum\limits_{l=0}^{n-k}A_{n-k-l}^{\alpha_{n}-1}\gamma_{l}$
$\displaystyle=$ $\displaystyle
A_{n-k}^{\alpha_{n}-1}\gamma_{0}+\sum\limits_{l=1}^{n-k}A_{n-k-l}^{\alpha_{n}-1}\gamma_{l}$
$\displaystyle=$ $\displaystyle
A_{n-k}^{\alpha_{n}-1}\left(1+\sum\limits_{l=1}^{n-k}\frac{A_{n-k-l}^{\alpha_{n}-1}}{A_{n-k}^{\alpha_{n}-1}}\gamma_{l}\right)$
$\displaystyle\geq$ $\displaystyle
A_{n-k}^{\alpha_{n}-1}\left(1+\sum\limits_{l=1}^{n-k}\frac{q_{n-k-l}}{q_{n-k}}\gamma_{l}\right)$
$\displaystyle=$
$\displaystyle\frac{A_{n-k}^{\alpha_{n}-1}}{q_{n-k}}\left(q_{n-k}+\sum\limits_{l=1}^{n-k}q_{n-k-l}\gamma_{l}\right)$
$\displaystyle=$
$\displaystyle\frac{A_{n-k}^{\alpha_{n}-1}}{q_{n-k}}\sum\limits_{l=0}^{n-k}q_{n-k-l}\gamma_{l}=0.$

Hence,

(17)
$\sum\limits_{k=0}^{n}\left|t_{k,n}\right|=\sum\limits_{k=0}^{n}t_{k,n}=1$

and we arrive at the condition $\left(c\right)$.

Now, let us choose the sequnce $(q_{n})$ as follows:

$q_{n}=\frac{1}{n+1}$

Then all the required conditions for $(q_{n})$ are fulfilled. Moreover, the
corresponding Nörlund means coincides with the Logarithmic means. Hence, from

(18)
$\sigma_{n}^{\alpha_{n}}(\mathbf{S})=\sum\limits_{k=0}^{n}t_{k,n}L_{k}\left(\mathbf{S}\right)$

we immediately arrive at the required assertion. ∎

The connection between divergence of the given sequence and the divergence of
the sequence $\sigma_{n}^{\alpha_{n}}(\mathbf{S})$ was established by
Tetunashvili [33]. In particularly, the following theorem was proved.

######  Theorem T.

Let $\left\\{\alpha_{n}\right\\}$ be a sequences such that

$0<\alpha_{n}\leq\frac{c}{\ln n},\quad\text{ where }0<c<\ln 2\text{ and }n>m.$

for some $m>0$. Then for any sequence
$\mathbf{S}:=\left\\{s_{n}:n\in\mathbb{N}\right\\}$:

$\sup\limits_{n}\left|s_{n}\right|=+\infty,$

implies

$\sup\limits_{n}\left|\sigma_{n}^{a_{n}}\right|=+\infty.$

From the previous Theorem 5 and Theorem T we obtain next result if one takes
$\alpha_{n}=\frac{1}{n+1}$.

######  Theorem 6.

Let $\mathbf{S}:=\left\\{s_{n}:n\in\mathbb{N}\right\\}$ be a sequence of real
numbers such that

$\sup\limits_{n}\left|s_{n}\right|=+\infty.$

Then

$\sup\limits_{n}\left|L_{n}\left(\mathbf{S}\right)\right|=+\infty.$

## Conflicts of Interest

The authors declare that they have no conflicts of interest.

## Declarations

### Ethical Approval

Not applicable

### Authorship Contribution

All authors contributed equally.

## References

  * [1] T. Akhobadze and Sh. Zviadadze.  A note on the generalized Cesáro means of trigonometric Fourier series.  Izv. Nats. Akad. Nauk Armenii Mat., 54(5):3–10, 2019.
  * [2] S. Bochkarev.  A Fourier series in an arbitrary bounded orthonormal system that diverges on a set of positive measure.  Matematicheskii Sbornik, 140(3):436–449, 1975.
  * [3] S. V. Bochkarev.  Everywhere divergent Fourier series in the Walsh system and in multiplicative systems.  Uspekhi Mat. Nauk, 59(1(355)):103–124, 2004.
  * [4] V. M. Bugadze.  Divergence of Fourier-Walsh series of bounded functions on sets of measure zero.  Mat. Sb., 185(7):119–127, 1994.
  * [5] V. V. Buzdalin.  Unboundedly diverging trigonometric Fourier series of continuous functions.  Mat. Zametki, 7:7–18, 1970.
  * [6] F. Di Plinio.  Lacunary fourier and walsh–fourier series near l^ 1 l 1.  Collectanea mathematica, 65:219–232, 2014.
  * [7] Y. Q. Do and M. Lacey.  On the convergence of lacunacy Walsh-Fourier series.  Bull. Lond. Math. Soc., 44(2):241–254, 2012.
  * [8] G. Gát.  Almost everywhere convergence of Fejér and logarithmic means of subsequences of partial sums of the Walsh-Fourier series of integrable functions.  J. Approx. Theory, 162(4):687–708, 2010.
  * [9] G. Gát.  Cesàro means of subsequences of partial sums of trigonometric Fourier series.  Constr. Approx., 49(1):59–101, 2019.
  * [10] G. Gát.  Almost everywhere divergence of cesáro means of subsequences of partial sums of trigonometric fourier series.  Mathematische Annalen, 2023.
  * [11] G. Gát and U. Goginava.  Uniform and $L$-convergence of logarithmic means of Walsh-Fourier series.  Acta Math. Sin. (Engl. Ser.), 22(2):497–506, 2006.
  * [12] G. Gát and U. Goginava.  On the divergence of Nörlund logarithmic means of Walsh-Fourier series.  Acta Math. Sin. (Engl. Ser.), 25(6):903–916, 2009.
  * [13] U. Goginava.  On the divergence of Walsh-Fejér means of bounded functions on sets of measure zero.  Acta Math. Hungar., 121(4):359–369, 2008.
  * [14] U. Goginava and G. Oniani.  On the divergence of subsequences of partial Walsh-Fourier sums.  Journal of Mathematical Analysis and Applications, 497(2):124900, 2021.
  * [15] R. Gosselin.  On the divergence of Fourier series.  In Proc. Amer. Math. Soc, volume 9, pages 278–282, 1958.
  * [16] G. H. Hardy.  Divergent series, volume 334.  American Mathematical Soc., 2000.
  * [17] J. Kahane and Y. Katznelson.  Sur les coefficients des séries de Fourier dont les sommes partielles sont positives sur un ensemble.  Studia Math., 44:555–562, 1972.
  * [18] I. B. Kaplan.  Cesàro means of variable order.  Izv. Vysš. Učebn. Zaved. Matematika, 1960(5 (18)):62–73, 1960.
  * [19] G. A. Karagulyan.  Divergence of general operators on sets of measure zero.  Colloq. Math., 121(1):113–119, 2010.
  * [20] Y. Katznelson.  Sur les ensembles de divergence des séries trigonométriques.  Studia Mathematica, 26(3):305–306, 1966.
  * [21] K. Kazaryan and P. Ul’yanov.  On a problem of alexits with respect to divergence of orthogonal fourier series.  In Dokl. Akad. Nauk SSSR, volume 316, pages 542–546, 1991.
  * [22] A. Kolmogoroff.  Une série de Fourier-Lebesgue divergente presque partout.  Fundamenta mathematicae, 4(1):324–328, 1923.
  * [23] A. Kolmogoroff.  Une série de Fourier-Lebesgue divergente partout.  CR Acad. ScL Paris, 183, 1926.
  * [24] S. Konyagin.  Almost everywhere convergence and divergence of Fourier series.  In International congress of mathematicians, volume 2, pages 1393–1403. Citeseer, 2006.
  * [25] S. V. Konyagin.  On a subsequence of Fourier-Walsh partial sums.  Mat. Zametki, 54(4):69–75, 158, 1993.
  * [26] S. V. Konyagin.  Divergence everywhere of subsequences of partial sums of trigonometric Fourier series.  Proc. Steklov Inst. Math., (Function Theory, suppl. 2):S167–S175, 2005.
  * [27] S. V. Konyagin.  Almost everywhere convergence and divergence of Fourier series.  In International Congress of Mathematicians. Vol. II, pages 1393–1403. Eur. Math. Soc., Zürich, 2006.
  * [28] S. V. Konyagin.  Almost everywhere divergence of lacunary subsequences of partial sums of Fourier series.  Proc. Steklov Inst. Math., 273(suppl. 1):S99–S106, 2011.
  * [29] V. Lie.  On the pointwise convergence of the sequence of partial fourier sums along lacunary subsequences.  Journal of Functional Analysis, 263(11):3391–3411, 2012.
  * [30] F. Móricz and A. H. Siddiqi.  Approximation by Nörlund means of Walsh-Fourier series.  J. Approx. Theory, 70(3):375–389, 1992.
  * [31] S. Stechkin.  On convergence and divergence of trigonometric series.  Uspekhi Mat. Nauk, 6(2):148–149, 1951.
  * [32] E. Stein.  On limits of sequences of operators.  Annals of Mathematics, pages 140–170, 1961.
  * [33] Sh. Tetunashvili.  On divergence of Fourier series by some methods of summability.  J. Funct. Spaces Appl., pages Art. ID 542607, 9, 2012.
  * [34] V. Totik.  On the divergence of fourier series.  Publ. Math. Debrecen, 29(3-4):251–264, 1982.
  * [35] A. Zygmund.  Trigonometric series, volume 2.  Cambridge university press, 2002.