Toeplitz package

[ivan-pi]

In regard to

I have not found any previous version of this algorithm and it is entirely my derivation. If I have reinvented the wheel, there being nothing new under the sun, it would be interesting to see the previous attribution.

There is package of routines on Netlib named TOEPLITZ: https://netlib.org/stoeplitz/

The TOEPLITZ package is a collection of Fortran subroutines for the numerical solution of systems of linear equations with coefficient matrices of Toeplitz or circulant form.

It was written in 1983 by a joint group of mathematicians and computer scientists from Argonne National Laboratory (USA) and the Moscow State University (USSR).

The full report can be found here: https://digital.library.unt.edu/ark:/67531/metadc283582/m1/1/

[ivan-pi]

Google Scholar finds only 28 works which cite the package: https://scholar.google.com/scholar?cites=9860393319970461535&as_sdt=2005&sciodt=0,5&hl=en&oi=gsb

[mostlyharmlessone]

hey ivan-pi, thanks for commenting. From what I understand, Toeplitz matrices have constant diagonals, with circulant matrices being a subset with all the diagonals based on a single vector. The matrices I define here are not Toeplitz matrices, as each aij is allowed to be unique.

[ivan-pi]

The report mentions four classes of matrices supported by the package:

• Toeplitz matrices (T-matrices)
• Circulant matrices (C-matrices)
• General matrices (G-matrices)
• Column-circulant matrices

Maybe yours fit under general?

I admit I haven't read the report or have a use for the routines. Your library appeared under my GitHub suggestions and I remembered seeing the package on Netlib with a similar scope.

A friendlier PDF version of the report can be found on the CERN server: https://cds.cern.ch/record/151420/files/CM-P00069203.pdf

[mostlyharmlessone]

I appreciate the pdf, particularly as the Google Scholar links are mostly behind pay walls. I read through the section on G-matrices, p.28 or so.

What the Toeplitz package calls a G-matrix is simply an NxN full matrix, not banded or sparse, and their routine simply provides an interface to LINPACK for those matrices. The matrices dealt with by my algorithm are general in the sense that the diagonals have no restrictions like Toeplitz (constant diagonals), circulant (square Toeplitz with identical but rotated diagonals) but are banded, ideally sparse, with NxM unique numbers in a NxN matrix.

Certainly if N=M then mine is a not very efficient solver for general square matrices: that is not the intent. On the other hand, the LINPACK routines referenced by the Toeplitz package for G-matrices take no advantage of the sparseness of the matrix. Note that LINPACK has been superseded by LAPACK.

The closest I could find were the LAPACK routines for general banded matrices, which I use as a comparison in my test suite, but those banded matrices are not cyclic; ie they have no non-zero elements outside the central band, whereas my matrices have the band "wrap-around" so that the off diagonal corners have non-zero elements.

[ivan-pi]

Right. Is there an application where such matrices originate from?

I can imagine something like a finite-difference method on a non-equidistant periodic grid would have tridiagonal cyclic matrix. In case of a bi-periodic grid it would be somehow block-band diagonal.

I think also a periodic cubic-spline with non-equidistant knots would have such a matrix. If I'm not mistaken there is an example in this book); it uses the Shermann-Morrison formula to deal with the wrap-around, but I don't know how generally applicable that is.

[mostlyharmlessone]

Well yes, as I point out in my notes, and as you point out, a finite difference equation with 3 points makes a triadiagonal matrix, and with periodic boundary conditions it is cyclic. A finite difference equation and periodic b.c. with 5 points makes a banded cyclic matrix with band width of 5. Irregular grids or non-constant or anisometric material properties, such as are common with finite element models, can lead to non-constant diagonals. In any case, it was an interesting problem, for me at least, regardless of applicability.

A periodic cubic spline with non-equidistant knots certainly leads to the tridiagonal version with non-constant diagonals. ("Numerical Recipes in C/Fortran/Etc." by Press et al. also solves that case with the Sherman-Morrison(-Woodbury) formula, just as in the book you reference.) However the Sherman-Morrison formula has to be applied for each additional non-zero element. In the case of the triadiagonal matrix that is only twice, once for each off diagonal term. For a 5-banded case (eg. a quintic spline with a 5 diagonal banded matrix) it would have to be applied 6 times, for a seven banded matrix, 12 times, etc. It is certainly an approach to the problem, however my algorithm is very different from that.

I did compare the speed of my algorithm with a triadiagonal version + Sherman-Morrison and found them nearly identical. I haven't tried encoding a Sherman-Morrison version applied many times automatically in combination with LAPACK's band solver, so I can't say which would win in the general case. Sounds like a fun project! I also can't speak to the stability, or accuracy, of many repeated S-M-W applications.