LU decomposition

Stream: general

[Yaël Dillies]

How much work is there towards LU factorization? @Jiale Miao, you did Gram-Schmidt and it's quite similar in spirit I believe.

[Eric Wieser]

There's no generalization from matrices to linear maps that's useful here, is there?

[Yaël Dillies]

Doesn't upper/lower triangularity somehow translate to something in terms of endomorphisms? Or does it not because it is not preserved under a change of basis?

[Yaël Dillies]

Also, we docs#matrix.block_triangular_matrix, but docs#matrix.det_of_lower_triangular is not stated in terms of it :thinking:

[Yaël Dillies]

I would generalize the b argument to block_triangular_matrix to o → α where has_lt α, so that we can state upper triangularity as well.

[Eric Wieser]

Block triangularity and triangularity probably aren't particularly easy to state with one definition

[Eric Wieser]

In the same way that docs#matrix.diagonal and docs#matrix.block_diagonal and docs#matrix.block_diagonal' are all separate

[Yaël Dillies]

I can define lower_triangular and upper_triangular in terms of my proposed generalized block_triangular, though. They are
block_triangular M id and block_triangular M to_dual respectively.

[Eric Wieser]

Your generalization sounds good either way (although maybe restrict it to a linear order)

[Eric Wieser]

We also seem to be missing the result that block_diagonal is block triangular

[Yaël Dillies]

Do we want lower_triangular and upper_triangular as standalone definitions?

[Eric Wieser]

Defined on fin n only?

[Yaël Dillies]

No, on a general order, without more complications.

[Eric Wieser]

Maybe as abbreviations

[Yaël Dillies]

Sure, that sounds fine.

[Yaël Dillies]

This gets me thinking. Is there a more general version of LU decomposition for block triangular matrices?

[Kevin Buzzard]

Regular LU is surely stronger than any statement you can make about decomposition into block lower and upper matrices. Or am I missing something? The general vocabulary for this is that lower and upper triangular matrices are called Borel subgroupsof GL(n), and block upper/lower triangular matrices are called parabolic subgroups. The general story is very well understood and works for reductive groups.

[Yaël Dillies]

What is "regular LU"? I only know what I was lectured last term.

[Kevin Buzzard]

Normal LU

[Kevin Buzzard]

What statement about LU for block triangular matrices isn't weaker than what you already know in the upper triangular case?

[Yaël Dillies]

For example, can't you get a decomposition of the form
$\begin{pmatrix} I_m & 0 \\ A & I_n \end{pmatrix}\begin{pmatrix} B & C \\ 0 & D \end{pmatrix}$

[Yaël Dillies]

You trade off some nonzero coefficients in U in exchange of zero coefficients in L under the diagonal.

[Mario Carneiro]

Isn't that just regular LU for matrices over the ring of matrices?

[Mario Carneiro]

we should definitely have the isomorphism between block matrices and matrices of matrices because a lot of theorems about block matrices follow from that equivalence

[Kevin Buzzard]

Yael I guess so, although this is LU plus one line.

[Kevin Buzzard]

If you write the decomposition as (lower unipotent) (diag) (Weyl group element) (upper unipotent) then I guess this beats everything

[Eric Wieser]

Matrices of matrices require docs#dmatrix for the interesting cases

[Eric Wieser]

And even then, you run into trouble with ⬝ vs *

[Yaël Dillies]

So what should the mathlib version be? I'm in the mood right now.

[Michael Stoll]

If you consider matrices of matrices, then all blocks are forced to have the same size, I think, but block matrices are more general.

[Yaël Dillies]

Eric Wieser said:

Matrices of matrices require docs#dmatrix for the interesting cases

:point_of_information:

[Kalle Kytölä]

Yaël Dillies said:

Doesn't upper/lower triangularity somehow translate to something in terms of endomorphisms? Or does it not because it is not preserved under a change of basis?

I think Kevin referred to something of this flavor, but not fully explicitly from what I saw: upper triangular matrices correspond exactly to those linear maps that leave invariant each of the subspaces in the complete flag associated with the standard basis. So to me it would not be a bad idea to define flags <https://en.wikipedia.org/wiki/Flag_(linear_algebra)> and at least relate triangularity of matrices to them (if not even define triangularity via this) as: the corresponding linear map leaves each of the subspaces of the flag invariant.

Not requiring the flag to be complete corresponds to block triangular matrices (if I understood correctly what that referred to).

But my first instinct is that the Borel subgroups of $GL_n$ should not have anything to do with how triangularity is actually implemented (although both Kevin and Wikipedia mention this), since we definitely want also non-invertible triangular matrices; not only the linear maps that stabilize a flag but also linear maps that merely leave all of its subspaces invariant. (Maybe the Lie algebra $\mathfrak{gl}_n$ of $GL_n$ is closer, but still does not really seem to help with the implementation.) (EDIT: Ah, perhaps Kevin meant the Lie algebra of the stabilizer subgroup of the flag! Anyways, I think those should be theorems, not a part of hte definition...)

And I think it would be great if one could handle matrices of matrices conveniently, but others have surely more insight into the feasibility of it.

[Yaël Dillies]

I am currently trying to generalize docs#matrix.det_of_block_triangular_matrix over at #14035 and it's a pain. This is THE central calculation, so I really don't know what to think of the matrix-free approach, especially because it won't be coordinate-free anyway.

[Kalle Kytölä]

A flag of subspaces is indeed specifying something, but not as much as coordinates/basis. And of course the complete flag corresponding to the standard basis is specifying a lot (although still not the normalizations of the basis vectors), but a special case naturally specifies more.

Anyways, I don't want to weigh in too much, just wanted to be more explicit about what was left a bit implicit before. Eventually there will be results saying that the matrix of a linear map that leaves the subspaces in the flag invariant is upper triangular in a basis that is adapted to the flag, and vice versa.

And I can easily believe the determinants are causing enough trouble without inviting more... :grinning_face_with_smiling_eyes:

[Kalle Kytölä]

Something that might become easier with the flag approach is, for example: the composition of upper triangular linear maps is upper triangular.

[Johan Commelin]

@Yaël Dillies already introduced flags a little while ago (in the more general context of ordered sets). So I think we can just talk about flag (submodule k V). Of course there is not linear-algebra-specific API for it yet.

[Yaël Dillies]

I think there is a slight disparity in the names. An order theory flag is a linear algebra complete flag, and we do not have linear algebra flags (aka partial flags yet). But it's super easy to add them.

[Eric Wieser]

docs#flag

[Eric Wieser]

Don't we only need _complete_ flags for triangular matrices?

[Yaël Dillies]

I think so, but not for block triangular matrices.

[Yaël Dillies]

I love emphasizing everything today.

[Yaël Dillies]

The relevant sentence on Wikipedia is

In the more abstract setting of incidence geometry, which is a set having a symmetric and reflexive relation called incidence defined on its elements, a flag is a set of elements that are mutually incident.[1] This level of abstraction generalizes both the polyhedral concept given above as well as the related flag concept from linear algebra.

[Yaël Dillies]

As I understand it, a partial flag is a chain that contains ⊥ and ⊤. What do you think of this as a definition?

[Damiano Testa]

Why do you want to automatically include top and bottom?

[Yaël Dillies]

This is what Wikipedia does, but I agree it's a bit redundant because we can always formally add them to any chain when they exist and it will still be a chain.

[Damiano Testa]

(I genuinely do not know: in "normal" maths it is so easy to say "oh, yes, ignore/include the zero vector space".)

[Yaël Dillies]

Not including the top/bottom element by default has a priori only one consequence, which is that "adding them later" is not injective (because the chain might already contain them). I don't know if this is a bug or a feature.

[Yaël Dillies]

But I like the idea that linear algebra flags are just chain (submodule 𝕜 V).

[Yaël Dillies]

Oh shoot, I didn't actually add chain. Let me do that quick.

[Damiano Testa]

Chains are very prominent:

• chains of ideals --> Krull dimension,
• chains of non-zero-divisors --> regularity,
• chains of normal subgroups --> Jordan-Hölder filtrations,
• chains of normal subgroups with abelian quotients --> solvability,
• ...

[Kevin Buzzard]

Don't chains have < not <=? If so then you can't add bot randomly to a chain

[Yaël Dillies]

The chains I'm talking about are sets, not indexed families. Indexed chains are just strictly monotone functions (docs#strict_mono).

[Kalle Kytölä]

Damiano Testa said:

Why do you want to automatically include top and bottom?

One reason might be the following. If $V_0=\{0\} \subsetneq V_1 \subsetneq \cdots \subsetneq V_n = V$, and a linear map $L \colon V \to V$ respects this flag (i.e. $L V_j \subset V_j$ for each $j$), then $L$ induces also endomorphisms in each of the successive quotients $V_j / V_{j-1} \to V_j / V_{j-1}$ --- these are essentially the "diagonal blocks". It seems very natural to include the first quotient $V_1 / V_0 = V_1 / \{0\} \cong V_1$ and the last one, $V_n/V_{n-1} = V / V_{n-1}$, so as to try to capture all information about the diagonal blocks of $L$.

For example a linear map is strictly upper triangular w.r.t. this flag if it respects the flag (each $V_j$ is an invariant subspace) and the mappings it induces $V_j / V_{j-1} \to V_j / V_{j-1}$ are all zero.

[Damiano Testa]

Sure, I do not object that complete flags, or also simply flags that contain the top and bottom element are special and important. I was wondering whether it was possibly useful to not require this as part of the definition (and maybe making a further definition of a "strict" flag as one containing top and bottom).

[Damiano Testa]

For instance, I can see wanting to consider a chain of prime ideals in a commutative ring over a field to be also a chain of vector spaces. Except that the whole space is not a prime ideal and the zero vector space need not be prime, if the ring contains non-zero zero-divisors.

I think that it might make sense to not force top and bottom to be parts of flags in general.

[Yaël Dillies]

I need something cleared out, however. Do you want flags to be sets or indexed families?

[Damiano Testa]

I think that most of the times, I've looked at finite flags and was often interested in taking successive quotients. This makes me think that indexed families might be more appropriate. However, again, the conversion is not such a big deal in maths, so I may be wrong about this.

[Damiano Testa]

In fact, I have to think hard about a situation where a "flag" was not simply a finite sequence of inclusion of subobjects, where the index set was 0, 1,..., n.

[Yaël Dillies]

Okay, so indexed families. In that case, docs#flag won't be terribly useful.

[Jireh Loreaux]

@Damiano Testa you haven't been doing enough infinite dimensional operator theory :smiley:. There's a whole theory of nest algebras, which are basically (not selfadjoint) operator algebras with a flag of invariant subspaces. However, the flag need not be finite, nor even countable, nor even well-ordered. There are lots of reasons to not be too restrictive about the index type of a flag.

[Damiano Testa]

To be fair, I have seen infinite flags: I love polynomials and I do not usually give a bound to their degree! So, I will accept nat as an indexing set. :stuck_out_tongue_wink:

[Yaël Dillies]

But are you flags still indexed, Jireh?

[Jireh Loreaux]

Sometimes it is useful to have an index, yes. Another classic example is a resolution of the identity corresponding to a selfadjoint operator, where you want the index set to be ℝ, or some interval therein.

[Yaël Dillies]

So basically a flag is a strictly monotone submodule-valued function?

[Jireh Loreaux]

yes, that's how I'm familiar with the term.

[Damiano Testa]

Me too, although the source of the function is a very tame linearly ordered set for me! :smile:

[Yaël Dillies]

Something like this?

structure partial_flag (α 𝕜 E : Type*) [preorder α] [semiring 𝕜] [add_comm_monoid E] [module 𝕜 E] :=
(to_fun : α → submodule 𝕜 E)
(strict_mono' : strict_mono to_fun)

[Yaël Dillies]

I don't know whether it warrants a new definition. Do we ever consider the type of flags?

[Jireh Loreaux]

I wouldn't define it now unless someone needs it (I don't currently). Otherwise we'll probably end up needing to redesign various bits of it when the time comes.

[Damiano Testa]

In algebraic geometry, the "set" of all flags carries the structure of an algebraic variety. The set of maximal flags is called the flag variety. These are very important, but probably to get the design right, there is the need to see where the algebraic geometry library ends up!

[Yaël Dillies]

Oh, flag might be the correct way to model the flag variety, because it's less nice to express that an indexed family doesn't "skip" elements (you can say that the range is order connected, I suppose).

[Damiano Testa]

If you look at flags that are not necessarily maximal, in the context of the flag variety, you get a (possibly) disconnected space, whose connected components are the ones where you keep track of which subset of {0,1,...,dim V} is the set of dimensions of elements of your flag.

[Eric Wieser]

Do we have the statement about docs#flag of set_range f where f is strict_mono?