Implement left_orth/right_orth

[mtfishman]

To do:

• Implement block sparse left_polar/right_polar.
• Add tests.

Future work:

• Implement block sparse lq_compact/lq_full.

Discussion points:

• Move out-of-place implementation of left_orth and right_orth to MatrixAlgebraKit.jl?

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 78.26%. Comparing base (13d0332) to head (be112c8).
Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #122      +/-   ##
==========================================
+ Coverage   77.61%   78.26%   +0.65%     
==========================================
  Files          32       33       +1     
  Lines        1604     1652      +48     
==========================================
+ Hits         1245     1293      +48     
  Misses        359      359              

Flag	Coverage Δ
docs	10.43% <0.00%> (-0.32%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[mtfishman]

@lkdvos going into this, I was interested to see how much of the logic in MatrixAlgebraKit.jl could be reused here. The primary issue I found was that left_orth! and right_orth! in MatrixAlgebraKit.jl are written in a way where they assume the output is pre-allocated:

• https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/blob/v0.2.0/src/interface/orthnull.jl#L72-L74
• https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/blob/v0.2.0/src/implementations/orthnull.jl#L54-L67

That isn't so straightforward for block sparse matrices and graded arrays since the block sizes and sectors may not be consistent across SVD, polar, and QR/LQ. That meant I had to repeat very similar logic in this PR to circumvent that. This PR makes me think that the MatrixAlgebraKit.jl logic is a bit too opinionated and could match the logic in this PR, I'm curious to hear what you think about that. Maybe MatrixAlgebraKit.jl could have a special codepath for StridedMatrix that allows pre-allocating the output.

[lkdvos]

Overall definitely looks reasonable to me.

As a global comment for the polar decomposition, the question is if we want to implement this as a blockwise polar, or by reimplementing the logic for computing the polar by doing a total svd and then reconstructing the polar decomposition.
Trying to think this through, I'm not sure which solution is better.
On the one hand, blockwise polar decomposition would require again to repeat the logic of constructing the correct axes, but it does seem more in the spirit of the blockwise algorithms.
On the other hand, your current implementation is quite simple and convenient, and probably the efficiency is more or less equal.
Just wanted to hear if you had more thoughts on this?

[lkdvos]

Is this not caught by the MatrixAlgebraKit implementation? If not, it probably should?

[mtfishman]

It's a bit subtle, the MatrixAlgebraKit.jl version is designed around the output being pre-allocated, so it has an extra argument for the output. That is tricky to support for block sparse matrices since the different backends (SVD, QR/LQ, polar) might have different block structures/sectors so it isn't straightforward to pre-allocate the output in that way.

We can move this code over to MatrixAlgebraKit.jl, that was the main thing I wanted your feedback on.

[lkdvos]

Ah I see now! This was indeed something that I at some point saw in MatrixAlgebraKit but then kind of forgot about in the end.
My initial reaction would be that I do think this should indeed be moved over there, such that the general codepath follows that of the @algdef-defined functions, with a little extra implementations that make it hard to simply use @algdef.

[mtfishman]

Ok, in that case I'll make a PR to MatrixAlgebraKit.jl, it would be nice to not have this logic here.

[mtfishman]

As a global comment for the polar decomposition, the question is if we want to implement this as a blockwise polar, or by reimplementing the logic for computing the polar by doing a total svd and then reconstructing the polar decomposition. Trying to think this through, I'm not sure which solution is better. On the one hand, blockwise polar decomposition would require again to repeat the logic of constructing the correct axes, but it does seem more in the spirit of the blockwise algorithms. On the other hand, your current implementation is quite simple and convenient, and probably the efficiency is more or less equal. Just wanted to hear if you had more thoughts on this?

It's a good question. I basically implemented it the current way because it is simpler and leverages the blockwise SVD. I think the leading scaling (in the limit of large blocks) should be the same so I'd lean towards keeping this version for the sake of code simplicity. Maybe we can switch over to the blockwise polar version when we change the logic of all of the blockwise factorizations to be based on first converting to a block diagonal matrix, in which case blockwise factorizations will be simpler to implement.

[mtfishman]

@lkdvos in the latest, I removed the implementation of polar since I moved it out to #125.

Also, inspired by our discussion above and on Slack, I simplified the logic here by overloading initialize_output for block sparse left_orth!/right_orth! to return nothing, so then that can get passed through the existing code logic in MatrixAlgebraKit. Then, I only overload the left_orth_f!/right_orth_f! functions. It feels slightly hacky but I think it is a reasonable way to go.

[lkdvos]

I definitely agree, the only thing I am wondering about is what happens if you actually do want to provide the output because you already have it. Do we just postpone handling that for now, and discard the provided outputs?

[mtfishman]

I definitely agree, the only thing I am wondering about is what happens if you actually do want to provide the output because you already have it. Do we just postpone handling that for now, and discard the provided outputs?

Yeah, for now it will just discard the provided outputs, which seems ok.

[mtfishman]

Actually, maybe I'll just make that case error for now so we know if we accidentally use that.

[mtfishman]

This is good to go from my end, any more comments?

[lkdvos]

I'm not entirely convinced by making that error - if you write generic code that is supposed to handle both dense and blocksparse arrays you would also not be able to provide the output. Realistically it doesn't matter that much yet though, so I'm okay with merging as is

[mtfishman]

I'm not entirely convinced by making that error - if you write generic code that is supposed to handle both dense and blocksparse arrays you would also not be able to provide the output. Realistically it doesn't matter that much yet though, so I'm okay with merging as is

I was thinking that too, but also the output isn't easy to construct for block sparse arrays anyway, so how easily can you write generic code across both cases? We can always remove the error if that case comes up, I just figure the error message will make it clear to us when that happens (and generic code can specify the output as nothing in the block sparse case anyway).