Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

block tridiagonal matrices #21

Closed

Conversation

maximerischard
Copy link

This is a quick attempt at implementing block tridiagonal matrices. I'm working my way towards symmetric block tridiagonal matrices, and their Cholesky decomposition (which would be a block lower bi-diagonal matrix? not sure what to call it). The code is modeled after blockarray.jl, but restricted to two dimensions (I'm not sure what a tridiagonal matrix would look like in higher dimensions). Let me know if this is something you would be interested in including to BlockArrays.jl, and I can add some tests.

@maximerischard
Copy link
Author

@tkelman I'd be interested to hear if this is potentially useful to you, or if you've already implemented some of this, as I gather from (JuliaLang/LinearAlgebra.jl#136) that you've dealt with this kind of object before.

@codecov-io
Copy link

codecov-io commented Jul 10, 2017

Codecov Report

Merging #21 into master will decrease coverage by 17.86%.
The diff coverage is 0%.

Impacted file tree graph

@@             Coverage Diff             @@
##           master      JuliaLang/julia#21       +/-   ##
===========================================
- Coverage   82.98%   65.12%   -17.87%     
===========================================
  Files           6        7        +1     
  Lines         288      367       +79     
===========================================
  Hits          239      239               
- Misses         49      128       +79
Impacted Files Coverage Δ
src/blocktridiag.jl 0% <0%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4ddb504...f0ce697. Read the comment docs.

@KristofferC
Copy link
Collaborator

Thanks for the PR! I am traveling right now so will take a bit before I look at this. I would be interested if you have any comments regarding the API / code structure and where you think things could improve.

@maximerischard
Copy link
Author

Hi Kristoffer,
It seemed obvious that a lot of care went into designing the abstract types and data structures. That made it easy to implement the block tridiagonal without touching any of the original code. Thank you for building this package.

There's a few comments I would make:

  1. I'm not sure I like that the provided constructors return a block matrix with unallocated blocks. I would advocate for actually allocating the blocks by default, and letting the user copy! into them. It just feels a little fragile, and I've shot myself in the foot a few times already (resulting in segfaults). It could be possible to start with unallocated blocks, but it shouldn't be the default.
  2. It also seems like there should be checks somewhere that the blocks have the correct size (according to the BlockSizes object). If not in an inner constructor, at least in some of the outer constructors.
  3. The functions that use macros are hard to read, and they mostly all do the same thing: iterate over blocks. It would be nice to put that logic in just one place, so those functions are easier to understand. Maybe an iterator of block indices? Not sure exactly...

for example `Matrix{Float64}.
"""
struct BlockTridiagMatrix{T, R <: AbstractMatrix{T}} <: AbstractBlockMatrix{T}
diagl::Vector{R}
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

An alternative would be to store the blocks in a BandedMatrix from the BandedMatrices.jl package, which would then take care of the logic at the block level. But then I don't know if introducing a dependency on BandedMatrices.jl is a good idea.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if the name should be TridiagBlockMatrix?

@maximerischard
Copy link
Author

I just saw that there's a fairly mature implementation of banded block matrices in ApproxFun.jl, and a placeholder for it in JuliaMatrices. @dlfivefifty, do you have any thoughts on where banded block matrices should belong? It seems like there's a lot of common functionality in BlockArrays.jl and ApproxFun.jl that would ideally be deduplicated.

@dlfivefifty
Copy link
Member

See discussion in #22
#23
I agree that the two almost-identical approaches should be merged. But it's not clear how to do some things I want to in the current BlockArrays.jl framework.

Note that while ApproxFun.BlockBandedMatrix Is relatively full featured, it's a bit hacky in places so will need to be fleshed out.

Do you think BlockTriadiagonal is distinctive enough that it should be its own type, a la SymTridiagonal? I think block tridiagonal matrices are especially important because of multivariate orthogonal polynomials.

@maximerischard
Copy link
Author

Thank you, I hadn't seen those issues. I'll go back and have a look through them. And no, I don't see any immediate reason why tridiagonal shouldn't be a special case of banded.

@KristofferC
Copy link
Collaborator

Note that I don't think this package has many users so if there are more convenient APIs we can simply switch. The stuff in this package is just what I happened to come up with for my use case.

@KristofferC
Copy link
Collaborator

Regarding point 3, the difficulty is to support arbitrary rank arrays. You could either do it with codegenning nested loops like I do with Cartesian macros or use some sort of recursion.

for example `Matrix{Float64}.
"""
struct BlockTridiagMatrix{T, R <: AbstractMatrix{T}} <: AbstractBlockMatrix{T}
diagl::Vector{R}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if the name should be TridiagBlockMatrix?


# Auxilary outer constructors
function BlockTridiagMatrix{T, R <: AbstractArray{T}
}(diagl::Vector{R},
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In my opinion, these linebreaks before the } looks a bit weird and I haven't seen Julia code using it before.


function BlockTridiagMatrix{T, R <: AbstractMatrix{T}}(::Type{R}, block_sizes::BlockSizes{2})
n_blocks = nblocks(block_sizes)
n_blocks[1] == n_blocks[2] || throw("expect same number of blocks in both dimensions")
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should probably throw an ArgumentError. I also prefer to write erros in "must"-form e.g the number of blocks must be the same in both dimensions.

# otherwise return a freshly-baked
# matrix of zeros (with a warning
# because that's dumb)
warn(@sprintf("""The (%d,%d) block of a block tridiagonal matrix
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't like warnings like this. It should either be an error or be allowed without warning. Potential wastefulness should go in the documentation of the type. You don't get a warning if you access a zero element of a sparse matrix.

return arr
end

function Base.copy!{T, R<:AbstractArray{T}, M<:AbstractArray{T}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Doesn't need to be done here but there should be a way where the code here doesnt have to be so similar to the one working for BlockArray.

@dlfivefifty
Copy link
Member

@maximerischard If you are still interested in this pull request, it makes more sense to live in the (currently in-development) repository https://github.com/JuliaMatrices/BlockBandedMatrices.jl

This way a common interface for talking about block-banded matrices will be used.

I'll close this pull request now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants