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.

Sign up for GitHub

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

Jump to bottom

Constructors from Kronecker Products? #62

Open
ChrisRackauckas opened this issue Feb 9, 2020 · 9 comments
Open

Constructors from Kronecker Products? #62

ChrisRackauckas opened this issue Feb 9, 2020 · 9 comments

Comments

@ChrisRackauckas
Copy link

Usually these types of matrices from from PDEs where they can be constructed from the Kronecker product of the actions along each dimension, like:

using SparseArrays
Iy = SparseMatrixCSC(I,N,N)
Ix = SparseMatrixCSC(I,N,N)
fJ = ones(3,3)
Dz = [1 0 0
      0 0 0
      0 0 0]
A = kron(Dz,Iy,sparse(Mx)) + kron(Dz,sparse(My),Ix) + kron(fJ,Iy,Ix)

Could there be a way like that last statement to directly form a BlockBandedMatrix or BandedBlockBandedMatrix?
@dlfivefifty
Copy link
Member

What's Mx and My?
Note that there is support for Kronecker products via LazyBandedMatrices:

using LazyBandedMatrices, BandedMatrices, BlockBandedMatrices, FillArrays, LazyArrays
N = 2
Ix = Iy = BandedMatrix(Eye(N,N))
julia> BandedBlockBandedMatrix(Kron(Ix,Iy))
2×2-blocked 4×4 BandedBlockBandedMatrix{Float64,BlockArrays.PseudoBlockArray{Float64,2,Array{Float64,2},Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}},BlockArrays.BlockedUnitRange{Array{Int64,1}}}:
 1.0   ⋅   │   ⋅    ⋅ 
  ⋅   1.0  │   ⋅    ⋅ 
 ──────────┼──────────
  ⋅    ⋅   │  1.0   ⋅ 
  ⋅    ⋅   │   ⋅   1.0

@ChrisRackauckas
Copy link
Author 

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

@ChrisRackauckas
Copy link
Author

Seems like there's some missing methods:

using LinearAlgebra, LazyBandedMatrices, BandedMatrices, BlockBandedMatrices, FillArrays, LazyArrays

N = 32
Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

Ix = BandedMatrix(Eye(N,N))
Iy = BandedMatrix(Eye(N,N))
Ix = SparseMatrixCSC(I,N,N)
fJ = ones(3,3)
Dz = [1 0 0
      0 0 0
      0 0 0]
A = BlockBandedMatrix(Kron(Dz,Iy,sparse(Mx))) +
    BlockBandedMatrix(Kron(Dz,sparse(My),Ix)) +
    BlockBandedMatrix(Kron(fJ,Iy,Ix))
MethodError: no method matching LazyBandedMatrices.BlockKron{Float64,A,B} where B where A(::Array{Int64,2}, ::BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}, ::SparseMatrixCSC{Float64,Int64})
Closest candidates are:
  LazyBandedMatrices.BlockKron{Float64,A,B} where B where A(::AA, ::BB) where {T, AA, BB} at C:\Users\accou\.julia\packages\LazyBandedMatrices\l1EHv\src\LazyBandedMatrices.jl:331
LazyBandedMatrices.BlockKron(::ApplyArray{Float64,2,typeof(kron),Tuple{Array{Int64,2},BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}},SparseMatrixCSC{Float64,Int64}}}) at LazyBandedMatrices.jl:333
BlockSkylineMatrix{T,Array{T,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockArrays.BlockedUnitRange{Array{Int64,1}},BlockArrays.BlockedUnitRange{Array{Int64,1}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},BandedMatrix{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}} where T(::ApplyArray{Float64,2,typeof(kron),Tuple{Array{Int64,2},BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}},SparseMatrixCSC{Float64,Int64}}}) at LazyBandedMatrices.jl:359
top-level scope at test.jl:102

@ChrisRackauckas
Copy link
Author

This example is trying to do SciML/SciMLBenchmarks.jl#91 but create a BandedBlockBandedMatrix, instead of doing a flat sparse matrix from SparsityDetection.jl

@dlfivefifty
Copy link
Member

You don’t want dense arguments (e.g. Dz ) in the construction of a BandedBlockBandedMatrix

Also, is a triple kron correspondent to a 3D PDE? That would require a more sophisticated data structure...

@ChrisRackauckas
Copy link
Author

It's a system of 3 2D PDEs.

@ChrisRackauckas
Copy link
Author

It's essentially

u_t = u_xx + u_yy + f1(u,v,w)`
v_t = f2(u,v,w)
w_t = f3(u,v,w)

The dense matrix is then the sparsity pattern of f.

@dlfivefifty
Copy link
Member

Yes that still requires a more sophisticated data structure: I believe if we interlace the rows one gets a banded-block-banded structure, which is what ApproxFun does but with a hacky implementation (InterlaceOperator). I had started a while ago to make an Interlace matrix but never wrapped it up

@ChrisRackauckas
Copy link
Author

ChrisRackauckas commented Feb 9, 2020 •
edited
Loading

Got it. I thought it ended up as a block banded matrix in the end (but not banded block banded).

It should be blocks of the block banded matrix of the 2D Laplacian operator.

@ChrisRackauckas ChrisRackauckas mentioned this issue Feb 26, 2020
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@dlfivefifty @ChrisRackauckas