Create blockcholesky.jl #126

felixw98 · 2020-07-14T00:23:53Z

No description provided.

codecov · 2020-07-14T00:35:56Z

Codecov Report

Merging #126 (99f54d5) into master (0551a95) will increase coverage by 0.31%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #126      +/-   ##
==========================================
+ Coverage   92.29%   92.61%   +0.31%     
==========================================
  Files          14       15       +1     
  Lines        1376     1435      +59     
==========================================
+ Hits         1270     1329      +59     
  Misses        106      106

Impacted Files	Coverage Δ
src/blockcholesky.jl	`100.00% <100.00%> (ø)`
src/blocklinalg.jl	`96.53% <100.00%> (ø)`
src/blockproduct.jl	`93.10% <100.00%> (+0.51%)`	⬆️

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 0551a95...99f54d5. Read the comment docs.

dlfivefifty

Please add unit tests by creating a file test/test_blockcholesky.jl

dlfivefifty · 2020-07-14T09:39:42Z

src/blockcholesky.jl

+# Function for cholesky decomposition on block matries
+# 'cholesky' is the build-in one in LAPACK
+
+function Bcholesky(A::Symmetric{<:Any,<:BlockArray})


In src/BlockArrays.jl, add

import LinearAlgebra: cholesky, cholesky!

And then you can overload cholesky and cholesky!

dlfivefifty · 2020-07-14T09:47:57Z

src/blockcholesky.jl

+    chol_U = BlockArray{Float64}(zeros(size(A)), fill(size(A)[1]÷blocksize(A)[1], blocksize(A)[1]), fill(size(A)[1]÷blocksize(A)[1], blocksize(A)[1]))
+
+    # Initializing the first role of blocks
+    chol_U[Block(1,1)] = cholesky!(A[Block(1,1)]).U


When we change this to be properly in-place, we will want to replace this with cholesky!(Symmetric(getblock(parent(A),1,1))) which will modify A.

dlfivefifty · 2020-07-14T09:49:44Z

src/blockcholesky.jl

+    # Initializing the first role of blocks
+    chol_U[Block(1,1)] = cholesky!(A[Block(1,1)]).U
+    for j = 2:blocksize(A)[1]
+        chol_U[Block(1,j)] = chol_U[Block(1,1)]' \ A[Block(1,j)]


When this is in-place it will look something like:

P = parent(A) ldiv!(UpperTriangular(getblock(P,1,1))', getblock(P,1,j))

dlfivefifty · 2020-07-14T09:52:45Z

src/blockcholesky.jl

+    for i = 2:blocksize(A)[1]
+        for j = i:blocksize(A)[1]
+            if j == i
+                chol_U[Block(i,j)] = cholesky!(A[Block(i,j)] - sum(chol_U[Block(k,j)]'chol_U[Block(k,j)] for k in 1:j-1)).U


Here we will want to use mul! to do this in-place. Something like:

Pij = getblock(P,i,j) for k = 1:j-1 mul!(Pij,getblock(P,k,j)',getblock(P,k,j),-1.0,1.0) end cholesky!(Pij)

dlfivefifty · 2020-07-18T21:00:01Z

src/blockcholesky.jl

+
+"""
+Instructions and examples
+"""


Why are these strings included in the file? It's better to just add them as comments in the PR.

dlfivefifty · 2020-07-18T21:01:41Z

src/blockcholesky.jl

+1.changed all computations into in-place but there are still many allocations 
+which is highly depended on the number of blocks.
+For example, I took a block matrix A with n=5 and d=3 and it had 59 allocations. The build-in only have 5 instead.
+I will focus on reducing these allocations.


Note Julia v1.4 allocates for things that shouldn't, for example subviews, so you won't ever get 0 allocations.

This has been improved in Julia v1.5, but expecting 0 allocations is unrealistic at this point.

OK. I have played around with these today, finding that my function uses less memory but slower. Is this a trade-off or a result of inefficient algorithm I adopted?

dlfivefifty · 2020-07-18T21:03:04Z

I think the most important thing to add right now are unit tests.

felixw98 · 2020-07-18T23:57:05Z

I think the most important thing to add right now are unit tests.

Is that the test with build-in cholesky? I did that already on various block matrices and all passed.
These are the codes
A = BlockArray{Float64}(randn(100,100)+107I, fill(10,10), fill(10,10)); A = Symmetric(A)
MT = Matrix(A)

@test cholesky(A) ≈ cholesky(MT).U

dlfivefifty · 2020-07-19T05:29:44Z

I mean add a new file test/test_cholesky.jl with the unit tests.

dlfivefifty · 2020-07-20T07:11:36Z

Make sure you include the new test file in test/run tests.jl

dlfivefifty · 2020-07-21T07:13:49Z

Unfortunately the 5-argument mul! is not available in Julia v1.1. The easiest solution is to change it to ArrayLayouts muladd!:
mul!(C, A, B, α, b) becomes muladd!(α, A, B, β, C).

We also need it to return a Cholesky object. I think this is as simple as changing the return to Cholesky(chol_P, A.uplo, 0).

Note we should also support the case of Symmetric(A, :L), which would give a lower Cholesky factorisation.

dlfivefifty · 2020-07-21T16:05:29Z

Why did you close this?

felixw98 · 2020-07-21T16:13:02Z

Why did you close this?

Sorry, I clicked the wrong one by mistake. What you mean Symmetric(A, :L), if we want a lower Cholesky factorization, we can quote by .L? Or we shall have another analogue function for lower cases?

dlfivefifty · 2020-07-21T16:15:01Z

Yes. You would basically do the transpose of what we have, with Symmetric(getblock(...), :L) instead of Symmetric(getblock(...)).

See https://github.com/JuliaLang/julia/blob/07616a4e352de785db0a0a702c2a65fe905b864d/stdlib/LinearAlgebra/src/cholesky.jl#L217

dlfivefifty · 2020-07-23T07:29:20Z

src/blockcholesky.jl

        end
    end
+
+    return UpperTriangular(A)


Please have this return Cholesky(A, :U, 0). The 0 should actually be info which is determined from the cholesky! call, so probably what we want is

info = 0 for ... ... info = cholesky!(Symmetric(Pii)) iszero(info) || break ... end Cholesky(A, :U, info)

though we should double check what info actually means: that is, do we break and info indicates a problem, or do we continue as if it were SPD?

Similarly for lower below.

In the document of PosDefException(info), it says that the info field indicates the location of (one of) the eigenvalue(s) which is (are) less than/equal to 0. And the info is determined by the function LAPACK.potrf!(uplo, A). So we can do the check of SPD in a similar way with LinearAlgebra.cholesky

I think you mean indicates that there exists a negative eigenvalue, not that it finds one. Here is the doc for LAPACK.potrf!: http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html

INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the leading minor of order i is not positive definite, and the factorization could not be completed.

So I think we want

if !iszero(info) @assert info > 0 # should never have illegal values info += sum(blockaxes(A,2)[Block.(1:i-1)]) # add columns before we got to block to indicate correct minor break end

dlfivefifty · 2020-07-27T13:31:50Z

src/blockcholesky.jl


 cholesky(A::Symmetric{<:Real,<:BlockArray},
-    ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)
+    ::Val{false}=Val(false); check::Bool = false) = cholesky!(cholcopy(A); check = check)


 function b_chol!(A::BlockArray{T}, ::Type{UpperTriangular}) where T<:Real


Convention in Julia is that internal routines start with _. So perhaps rename this function _block_chol!?

dlfivefifty · 2020-07-27T13:32:33Z

src/blockcholesky.jl

-
-        # For the left blocks
-        for i = 2:n
+        for i = 1:n
            Pii = getblock(A,i,i) 
            for k = 1:i-1
                muladd!(-1.0, getblock(A,i,k), getblock(A,i,k)', 1.0, Pii)


1.0 should be one(T)

dlfivefifty · 2020-07-27T13:32:44Z

test/test_cholesky.jl

@@ -9,12 +9,20 @@ using BlockArrays, Test, LinearAlgebra
    B = BlockArray{Float64}(randn(55,55)+100I, 1:10, 1:10); B = Symmetric(B)


Add tests with Float32

dlfivefifty · 2020-08-25T12:46:43Z

src/blockbandedcholesky.jl

+
+
+
+function _blockbandedcholesky!(A::BlockArray{T}, ::Type{UpperTriangular}) where T<:Real


I don't understand why this is a separate function, not just block_chol! from beliw.

The only difference is that I added a check of trivial blocks from backwards. It actually do the same thing on dense matrices. We can have another parameter for activating the 'check'.

dlfivefifty · 2020-08-25T12:47:29Z

src/blockbandedcholesky.jl

+    @inbounds begin
+        for i = 1:n
+            Pii = getblock(A,i,i)
+            for k = 1:i-1


This should loop over colsupport(A.blocks, i) ∩ 1:i-1 I think

colsupport(_, A, j) = axes(A,1)

""""
colsupport(A, j)

gives an iterator containing the possible non-zero entries in the j-th column of A.
"""
colsupport(A, j) = colsupport(MemoryLayout(typeof(A)), A, j)

The code of colsupport(A, j) is just the axes(A,1), so this may not be useful.

It's overloaded for banded matrices: https://github.com/JuliaMatrices/BandedMatrices.jl/blob/3eb9c0d77f471db67b87b7dc5f9832c1d0beb0b2/src/generic/AbstractBandedMatrix.jl#L89

Since we only care about the blockbanded cases, maybe we can move the file into the package BlockBandedMatrices?
Then it will be convenient where all the prerequisite functions are included.

No, it should be here. It will do the right thing for both dense and banded formats.

Another thing is that why LAPACK.potrf!() do not work in-place on BlockBandedMatrices?
It gives the right output but the block will not be replaced by the decomposed matrix.

dlfivefifty · 2020-08-25T12:47:49Z

src/blockbandedcholesky.jl

+            end
+
+            k_start = n
+            while getblock(A,i,k_start) ≈ zeros(Float32, size(getblock(A,i,k_start))) && k_start > k_end


This is strange

felixw98 · 2020-08-31T20:28:44Z

julia> M
3×3-blocked 9×9 BlockSkylineMatrix{Float32,Array{Float32,1},BlockBandedMatrices.BlockSkylineSizes{Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},Fill{Int64,1,Tuple{Base.OneTo{Int64}}},BandedMatri
x{Int64,Array{Int64,2},Base.OneTo{Int64}},Array{Int64,1}}}:
 64.0  16.0   0.0  │  16.0   0.0   0.0  │    ⋅     ⋅     ⋅
 16.0  64.0  16.0  │   0.0  16.0   0.0  │    ⋅     ⋅     ⋅
  0.0  16.0  64.0  │   0.0   0.0  16.0  │    ⋅     ⋅     ⋅
 ──────────────────┼────────────────────┼──────────────────
 16.0   0.0   0.0  │  64.0  16.0   0.0  │  16.0   0.0   0.0
  0.0  16.0   0.0  │  16.0  64.0  16.0  │   0.0  16.0   0.0
  0.0   0.0  16.0  │   0.0  16.0  64.0  │   0.0   0.0  16.0
 ──────────────────┼────────────────────┼──────────────────
   ⋅     ⋅     ⋅   │  16.0   0.0   0.0  │  64.0  16.0   0.0
   ⋅     ⋅     ⋅   │   0.0  16.0   0.0  │  16.0  64.0  16.0
   ⋅     ⋅     ⋅   │   0.0   0.0  16.0  │   0.0  16.0  64.0

The block banded matrix like this and

julia> LAPACK.potrf!('U',getblock(M,1,1))
(Float32[8.0 2.0 0.0; 16.0 7.745967 2.065591; 0.0 16.0 7.7287345], 0)

But matrix M isn't changed. This won't happen on general blockarrays.

dlfivefifty · 2020-09-01T13:47:01Z

Can you try with view(M,Block(1,1))?

felixw98 · 2020-09-01T19:02:56Z

Can you try with view(M,Block(1,1))?

This works.
But in the following link, the rdiv!() takes the type of StridedMatrix{T} where it is StridedVecOrMat{T} in ldiv!().
Shouldn't these be the same?
https://github.com/JuliaLang/julia/blob/697e782ab86bfcdd7fd15550241fe162c51d9f98/stdlib/LinearAlgebra/src/triangular.jl#L766-L784

So far, rdiv!() cannot take view() as input since it isn't a subtype of StridedMatrix{T}.

dlfivefifty · 2020-09-02T08:20:41Z

Because you can't right-multiply a vector times a matrix

Just overload rdiv for BlockBandedBlock (defined here https://github.com/JuliaMatrices/BlockBandedMatrices.jl/blob/de57bfe664d62c6cbd354dd23d0c8b48a1aa8696/src/BlockSkylineMatrix.jl#L420)

But note I suggested you first try BlockArray{T,<:BandedMatrix} which is a different type than BlockBandedMatrix

felixw98 · 2020-09-03T14:31:22Z

BlockArray{T,<:BandedMatrix}

There isn't a type of BlockArray{T,<:BandedMatrix}, the only one I can get is BlockArray{T,2,<:BandedMatrix}.

julia> BlockArray{T,<:BandedMatrix} where T <: Real
ERROR: TypeVar in Vararg length must have bounds Union{} and Any
Stacktrace:
[1] top-level scope at REPL[81]:1

And for `BlockArray{T,2,<:BandedMatrix}`, `muladd!()` will have an error eg.
julia> A = BlockArray(A,fill(3,3),fill(3,3))
3×3-blocked 9×9 BlockArray{Float32,2,Array{BandedMatrix{Float32,Array{Float32,2},Base.OneTo{Int64}},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}:
 99.0099     0.595538     0.0       │    0.0        0.0         0.0       │   0.0         0.0         0.0
 -0.682299  99.6287      -0.254271  │    0.0        0.0         0.0       │   0.0         0.0         0.0
  0.0       -0.0365762  101.114     │    0.387621   0.0         0.0       │   0.0         0.0         0.0
 ───────────────────────────────────┼─────────────────────────────────────┼───────────────────────────────────
  0.0        0.0          0.887276  │  101.698     -1.04962     0.0       │   0.0         0.0         0.0
  0.0        0.0          0.0       │    0.335081  99.3336      1.1688    │   0.0         0.0         0.0
  0.0        0.0          0.0       │    0.0       -0.771624  100.44      │  -0.512122    0.0         0.0
 ───────────────────────────────────┼─────────────────────────────────────┼───────────────────────────────────
  0.0        0.0          0.0       │    0.0        0.0         0.944132  │  99.6787      0.122043    0.0
  0.0        0.0          0.0       │    0.0        0.0         0.0       │  -1.04835   101.11       -0.697174
  0.0        0.0          0.0       │    0.0        0.0         0.0       │   0.0        -0.331242  100.515

julia> ArrayLayouts.muladd!(-one(Float32), getblock(A,1,3)', getblock(A,1,3), one(Float32), getblock(A,3,3))
3×3 BandedMatrix{Float32,Array{Float32,2},Base.OneTo{Int64}}:
 0.0  0.0   ⋅
 0.0  0.0  0.0
  ⋅   0.0  0.0

while mul!() won't have this problem. The function replaces C with aAB instead of aAb + bC.

dlfivefifty · 2020-09-03T14:41:03Z

Meant BlockMatrix{T,<:BandedMatrix{<:Matrix}}. Your second example is BlockMatrix{T,<:Matrix{<:BandedMatrix}}}

felixw98 · 2020-09-04T01:28:13Z

Meant BlockMatrix{T,<:BandedMatrix{<:Matrix}}. Your second example is BlockMatrix{T,<:Matrix{<:BandedMatrix}}}

How to produce the first case, I have tried different ways but can only the later one.
And how to overload stride for BandedMatrix?

Create blockcholesky.jl

9c47f78

dlfivefifty requested changes Jul 14, 2020

View reviewed changes

dlfivefifty reviewed Jul 14, 2020

View reviewed changes

Add files via upload

030a054

dlfivefifty reviewed Jul 18, 2020

View reviewed changes

dlfivefifty and others added 2 commits July 19, 2020 07:14

overload cholesk!

b977d8d

Including the test file in REPL will do the test

0ff3e67

Including test_cholesky.jl in runtests.jl

dc54673

felixw98 closed this Jul 21, 2020

felixw98 reopened this Jul 21, 2020

Added lower Cholesky factorization

5605d8f

dlfivefifty reviewed Jul 23, 2020

View reviewed changes

Added PD check

bb7355b

dlfivefifty reviewed Jul 27, 2020

View reviewed changes

felixw98 added 2 commits July 27, 2020 18:09

Changed the function name and types

2d782fe

Added blockbanded case

72259da

dlfivefifty requested changes Aug 25, 2020

View reviewed changes

Rearrange

883da65

Capable on both BlockArray and BlockBandedArray

09a8a54

Changes

b6f7d3a

dlfivefifty mentioned this pull request May 4, 2022

Cholesky decompositon of spd BandedBlockBanded matrices JuliaLinearAlgebra/InfiniteLinearAlgebra.jl#104

Open

dlfivefifty added 9 commits June 14, 2022 15:05

Merge branch 'master' into pr/126

dc6e020

Update Choleksy

8793acd

tests pass

f994fc0

Support sparse block structure

44b1490

Update blockcholesky.jl

2d43bda

Fix tests on 1.8, use ArrayLayouts.ldiv!

025a4ce

Speed up BlockKron construction

04af833

Update test_blockproduct.jl

f1059ef

increase cov

99f54d5

dlfivefifty merged commit 32990fc into JuliaArrays:master Jun 16, 2022

timholy mentioned this pull request Jun 16, 2022

Support promotion against Base AbstractUnitRanges #213

Merged

		@@ -9,12 +9,20 @@ using BlockArrays, Test, LinearAlgebra
		B = BlockArray{Float64}(randn(55,55)+100I, 1:10, 1:10); B = Symmetric(B)




		function _blockbandedcholesky!(A::BlockArray{T}, ::Type{UpperTriangular}) where T<:Real

Create blockcholesky.jl #126

Create blockcholesky.jl #126

Conversation

felixw98 commented Jul 14, 2020

codecov bot commented Jul 14, 2020 • edited Loading

Codecov Report

dlfivefifty left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dlfivefifty commented Jul 18, 2020

felixw98 commented Jul 18, 2020

dlfivefifty commented Jul 19, 2020

dlfivefifty commented Jul 20, 2020

dlfivefifty commented Jul 21, 2020

dlfivefifty commented Jul 21, 2020

felixw98 commented Jul 21, 2020

dlfivefifty commented Jul 21, 2020

dlfivefifty Jul 23, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

felixw98 commented Aug 31, 2020

dlfivefifty commented Sep 1, 2020

felixw98 commented Sep 1, 2020 • edited Loading

dlfivefifty commented Sep 2, 2020

felixw98 commented Sep 3, 2020 • edited Loading

dlfivefifty commented Sep 3, 2020

felixw98 commented Sep 4, 2020

codecov bot commented Jul 14, 2020 •

edited

Loading

dlfivefifty Jul 23, 2020 •

edited

Loading

felixw98 commented Sep 1, 2020 •

edited

Loading

felixw98 commented Sep 3, 2020 •

edited

Loading