# Efficient Minors of a Toeplitz Matrix - Is There Demand?

**URL:** <https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318>\
**Category:** General Usage\
**Tags:** linearalgebra, contributing\
**Created:** [January 3, 2024, 11:07pm UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318 "2024-01-03T23:07:16Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [January 3, 2024, 11:07pm UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318/1 "2024-01-03T23:07:16Z")

</div>

I have had a use case there I need (the FFT of) the sequence of the primary minors of a square matrix `A`, e.g. schematically `fft([det(A[1:m, 1:m] for m in 1:size(A,1)])`. One would expect naively that this is `O(n^4)` in the size of `A`. In my case, `A` is a Toeplitz matrix. I found it a pleasant surprise that the Levison-Trench-Zohar algorithm (I found a nice summary [here](https://dspace.mit.edu/bitstream/handle/1721.1/4954/RLE-TR-538-20174000.pdf?sequence=1) ), not only computes a single minor in `O(n^2)`, it implicitly runs iteratively on minors, implying that a single run on `A` recovers the entire minor sequence I need. I have implemented this algorithm and am really pleased with the result.

My question is: would this be of interest for the ecosystem for me to contribute this implementation to, say, `ToeplitzMatrices.jl`? Specifically, the points I’d like to double-check is:

1. Is it worthwhile to specialize `LinearAlgebra.det` to use this? I infer that this algorithm might not necessarily stable (though positive definiteness is a sufficient condition), so this might not actually be desirable, or maybe a check-and-fallback would be ideal?
2. The Levinson algorithm is already implemented in `StatsBase`, and used by `ToeplitzMatrices`, so perhaps this leads to some duplication in the ecosystem? Yet that implementation is only for solving `A*x=b`, so I don’t see how I can use this implementation for determinants.
3. The real advantage came for my use case of needing the sequence of minors, which would not be exposed by just specializing `det`. Are there any thoughts on a proper interface to expose this idiomatically? One thought is to implement the UDL decomposition in the LTZ algorithm as a new factorization type (unless it already exists), and then read off my minors as the diagonal entries of D.
4. As a first-time contributor, if there is any SOP I should follow.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 4, 2024, 12:54am UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318/2 "2024-01-04T00:54:53Z")

</div>

IMO this probably shouldn’t exist as an optimization for `det`. `det` doesn’t have a lot of optimizations because it’s very rarely a thing you actually want to compute.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 4, 2024, 12:57am UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318/3 "2024-01-04T00:57:05Z")

</div>

> [@bremez](#):
>
> Is it worthwhile to specialize `LinearAlgebra.det` to use this?

It sounds like a potentially nice addition to ToeplitzMatrices.jl (though I’m not a primary maintainer of that package so my opinion is not dispositive).

> As a first-time contributor, if there is any SOP I should follow.

File a github pull request with your new function/method. Make sure it is documented and has a test.

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [January 10, 2024, 6:38pm UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318/4 "2024-01-10T18:38:49Z")

</div>

@stevengj For completion, is there a “standard convention” on what to do in case of ill-conditioned input? Silently return garbage? Throw `DomainError`? Under the hood switch to the general `det` implementation? (With `@info`?) I noticed things like `inv` of a singular matrix return an `LAPACKException` which I don’t find very informative.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2024, 7:02pm UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318/5 "2024-01-10T19:02:50Z")

</div>

> [@bremez](#):
>
> what to do in case of ill-conditioned input? Silently return garbage? Throw `DomainError`? Under the hood switch to the general `det` implementation?

What do you mean by ill conditioning in this context? An ill-conditioned matrix still has a perfectly well defined determinant. Of course, you can’t necessarily achieve a small relative error in this case (e.g. if the exact determinant is `0.0` and you return `1e-16`, that’s a relative error of infinity), so in that sense the determinant has a large relative condition number, but you should still be able to achieve a small error relative to other relevant scales like the norm^n of the matrix.

It’s also a matter of documentation and context. If there is a regime in which you know your algorithm gives inaccurate results, the ideal thing (besides improving your algorithm) is to switch to another algorithm that is more accurate in that case. But suppose that such an algorithm doesn’t exist, or is vastly more expensive (e.g. requires arbitrary precision arithmetic); in that case, depending on your expected users, you might silently return garbage but include a warning in the documentation, or throw an error (if you can identify when garbage results will occur), or do both (throw an error, but have a flag like `check=false` in the `lu` function that disables the error).

---

<div class="post-metadata">

**Author:** ![bremez](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bremez/32/38777_2.png) [@bremez](https://discourse.julialang.org/u/bremez)\
**Post date:** [January 10, 2024, 8:01pm UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318/6 "2024-01-10T20:01:20Z")

</div>

I meant a matrix that’s not positive-definite, which is a pre-requisite for the LTZ algorithm.

The algorithm implementation naturally can identify if this occurs: it successively computes the minors’ determinant, and then one would be zero (or test to be below a threshold). (The failure arises from a division by the latest minor determinant at each iteration.)

The flag option seems good, but differs from the semantics of `det`. It perhaps more strongly motivates an implementation of a UL factorization with LTZ, and then `det` of the factorization is trivial. [Side bar: are UL factorizations supported in Julia? My limited understanding is that the algorithm cannot be reformulated to produce an LU factorization.]

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 10, 2024, 8:27pm UTC](https://discourse.julialang.org/t/efficient-minors-of-a-toeplitz-matrix-is-there-demand/108318/7 "2024-01-10T20:27:34Z")

</div>

> [@bremez](#):
>
> Side bar: are UL factorizations supported in Julia?

Not directly, as far as I know, because they aren’t implemented in LAPACK. But you can get the UL factorization by doing [LU factorization of a permuted matrix](https://math.stackexchange.com/questions/1372166/a-ul-factorization).