Add elementary matrices constructors #21

rlschuller · 2021-09-12T06:26:59Z

As defined by Aitchison 1986.

A couple of observations regarding the choices of types:

Use Integer for i and j, following LinearAlgebra.jl conventions
Use Int for the variable inside the structure to avoid (re)allocation problems: it is used in F and G for the method getindex

The documentation was simplified to make it cleaner, but we can put the examples back.

juliohm

Thank you @rlschuller for the updated PR. It is very organized. I think we can do something more clever regarding the dimension of the matrices. Let's use I from LinearAlgebra as an example again:

julia> using LinearAlgebra

julia> I * [1,2,3]
3-element Vector{Int64}:
 1
 2
 3

It works even though we did not specify the dimension (=3). If you read their source you will see that they create a functor object that is dimensionless, and then when users type I(3) this creates a Diagonal instance. We can do something similar for our matrices so that when we use them in complicated expressions we don't need to type d everywhere.

Take a look at equations 4.26 -- 4.30 after the table in the book and imagine these expressions full of d, not very pleasant to read. We will use these equations to convert between the various types of covariance structure.

How to implement something similar to I without reinventing the wheel? Take a look at FillArrays.jl and BlockArrays.jl. For example I would export the functor J = JMatrix() and would write (J::JMatrix{T})(d) where {T} = Fill(one(T), d, d).

After we have this working with general dimension, we can check if the BlockArrays.jl approach introduces performance issues.

juliohm · 2021-09-12T11:55:46Z

src/CoDa.jl

@@ -26,6 +27,13 @@ export
  norm, dot, ⋅,
  distance,

+  # matrices
+  JColumn,
+  JMatrix,


Can you explain why you didn't export J and j as well?

Because of the conflicting styles - structures should start w/ upper case letters according to the manual.

There're is also the problem of j beeing a common indexation variable.

Do you have any suggestions for j alternatives? Maybe Jc?

Maybe for Julia this doesn't apply, but the variables I use most frequently in nested for loops are i, j and k. Wouldn't j for the column vector create confusion?

juliohm · 2021-09-12T12:03:49Z

test/matrices.jl

+  @test JColumn(d) == ones(d)
+  @test JMatrix(d) == ones(d, d)
+  @test F(d) == [I(d) -ones(d)]
+  @test maximum(abs.(G(D) - I(D)  + JMatrix(D)/ D)) < 1e-5


What about norm(G(D) - I(D) + J(D) / D, Inf) < 1e-5? The function norm is available in LinearAlgebra.

juliohm · 2021-09-12T12:04:00Z

test/runtests.jl

@@ -3,6 +3,7 @@ using DataDeps
 using RData
 using CSV
 using Test
+import LinearAlgebra: I


Suggested change

import LinearAlgebra: I

using LinearAlgebra

juliohm · 2021-09-12T12:05:51Z

test/matrices.jl

+  @test JColumn(d) == ones(d)
+  @test JMatrix(d) == ones(d, d)


I would rewrite these in terms of the shortcuts J and j after exporting these names.

juliohm · 2021-09-12T12:06:13Z

src/matrices.jl

+Base.IndexStyle(::Type{<:JColumn}) = IndexLinear()
+Base.getindex(A::JColumn{T}, i::Integer) where {T} = one(T)
+
+


Suggested change

juliohm · 2021-09-12T12:09:52Z

src/matrices.jl

+struct FMatrix{T} <: AbstractMatrix{T}
+  d::Int
+end
+FMatrix(d::Int) = FMatrix{Int}(d)


Notice that you are not consistent regarding the parameter type. For the Jmatrix you used Bool and here you used Int.

The elements in FMatrix can assume negative values.

The elements in GMatrix are fractional, so Float64 was chosen as the default. According to the Performance Tips section of the manual, the use of Real types should be avoided in containers, but for in my tests Real defaulted to Float64 and performed w/ the same speed.

Got it, thank you for clarifying it 💯

juliohm · 2021-09-12T12:10:18Z

src/matrices.jl

+end
+FMatrix(d::Int) = FMatrix{Int}(d)
+F(d::Int) = FMatrix(d)
+Base.size(A::FMatrix) = (A.d, A.d+1)


Suggested change

Base.size(A::FMatrix) = (A.d, A.d+1)

Base.size(F::FMatrix) = (F.d, F.d + 1)

I was unsure about how Julia would handle conflicts of names, if the arguments of the functions always overwrite the other variables w/ that name I agree that using F looks much more readable.

juliohm · 2021-09-12T12:12:08Z

src/matrices.jl

+F(d::Int) = FMatrix(d)
+Base.size(A::FMatrix) = (A.d, A.d+1)
+Base.IndexStyle(::Type{<:FMatrix}) = IndexCartesian()
+Base.getindex(A::FMatrix{T}, i::Integer, j::Integer) where {T} =  one(T)*((i==j) - (j==(A.d+1)))


Suggested change

Base.getindex(A::FMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*((i==j) - (j==(A.d+1)))

Base.getindex(F::FMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*((i == j) - (j == (F.d + 1)))

juliohm · 2021-09-12T12:13:32Z

src/matrices.jl

+end
+HMatrix(d::Int) = HMatrix{Int}(d)
+H(d::Int) = HMatrix(d)
+Base.size(A::HMatrix) = (A.d, A.d)


Suggested change

Base.size(A::HMatrix) = (A.d, A.d)

Base.size(H::HMatrix) = (H.d, H.d)

juliohm · 2021-09-12T12:14:44Z

src/matrices.jl

+H(d::Int) = HMatrix(d)
+Base.size(A::HMatrix) = (A.d, A.d)
+Base.IndexStyle(::Type{<:HMatrix}) = IndexCartesian()
+Base.getindex(A::HMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*(i==j) + one(T)


Suggested change

Base.getindex(A::HMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*(i==j) + one(T)

Base.getindex(::HMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*(i == j) + one(T)

rlschuller · 2021-09-12T13:33:10Z

I didn't realize that I was guessing the dimensions automatically, that's indeed a great feature.

juliohm · 2021-09-12T13:37:05Z

I didn't realize that I was guessing the dimensions automatically, that's indeed a great feature.

Yes, it is a nice side effect of multiple-dispatch in the language. The product A*x has a behavior that is determined by both arguments, in particular the second argument is what matters. So whenever I is passed as the first argument it is completely ignored. Similar expressions like x'*I*x can be implemented without specifying dimension of I. The beauty of Julia as I mentioned before :)

codecov-commenter · 2021-09-12T13:41:51Z

Codecov Report

Merging #21 (f0092ef) into master (54c3f44) will increase coverage by 2.15%.
The diff coverage is 100.00%.

❗ Current head f0092ef differs from pull request most recent head 4389bfa. Consider uploading reports for the commit 4389bfa to get more accurate results

@@            Coverage Diff             @@
##           master      #21      +/-   ##
==========================================
+ Coverage   82.10%   84.25%   +2.15%     
==========================================
  Files           3        4       +1     
  Lines          95      108      +13     
==========================================
+ Hits           78       91      +13     
  Misses         17       17

Impacted Files	Coverage Δ
src/matrices.jl	`100.00% <100.00%> (ø)`

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 54c3f44...4389bfa. Read the comment docs.

juliohm · 2021-09-12T15:11:47Z

I guess we could eliminate JColumn from the list temporarily. My impression is that most algebraic manipulations we will do can be expressed with just the matrices.

…

On Sun, Sep 12, 2021, 11:51 Rodrigo Schuller ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In src/CoDa.jl <#21 (comment)>: > @@ -26,6 +27,13 @@ export norm, dot, ⋅, distance, + # matrices + JColumn, + JMatrix, Maybe for Julia this doesn't apply, but the variables I use more frequently in nested for loops are i, j and k. Wouldn't j for the column vector create confusion? — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#21 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZQW3NZP2YTRQFSTZXQYPTUBS46DANCNFSM5D3YICZA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

As defined by Aitchison 1986.

juliohm · 2021-09-12T20:35:34Z

@rlschuller why you are closing the PRs instead of updating them in place? That way we loose all the suggestions in multiple places.

juliohm requested changes Sep 12, 2021

View reviewed changes

rlschuller force-pushed the rodrigo-task-2 branch 2 times, most recently from f0092ef to 6cb339d Compare September 12, 2021 20:32

Add elementary matrices constructors

4389bfa

As defined by Aitchison 1986.

rlschuller force-pushed the rodrigo-task-2 branch from 6cb339d to 4389bfa Compare September 12, 2021 20:33

rlschuller closed this Sep 12, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add elementary matrices constructors #21

Add elementary matrices constructors #21

rlschuller commented Sep 12, 2021

juliohm left a comment

juliohm Sep 12, 2021

rlschuller Sep 12, 2021

rlschuller Sep 12, 2021

rlschuller Sep 12, 2021 •

edited

Loading

juliohm Sep 12, 2021

juliohm Sep 12, 2021

juliohm Sep 12, 2021

juliohm Sep 12, 2021

juliohm Sep 12, 2021

rlschuller Sep 12, 2021

juliohm Sep 12, 2021

juliohm Sep 12, 2021

rlschuller Sep 12, 2021

juliohm Sep 12, 2021

juliohm Sep 12, 2021

juliohm Sep 12, 2021

rlschuller commented Sep 12, 2021

juliohm commented Sep 12, 2021 •

edited

Loading

codecov-commenter commented Sep 12, 2021 •

edited

Loading

juliohm commented Sep 12, 2021 via email

juliohm commented Sep 12, 2021

		Base.IndexStyle(::Type{<:JColumn}) = IndexLinear()
		Base.getindex(A::JColumn{T}, i::Integer) where {T} = one(T)

	Base.size(A::FMatrix) = (A.d, A.d+1)
	Base.size(F::FMatrix) = (F.d, F.d + 1)

	Base.getindex(A::FMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*((i==j) - (j==(A.d+1)))
	Base.getindex(F::FMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*((i == j) - (j == (F.d + 1)))

	Base.size(A::HMatrix) = (A.d, A.d)
	Base.size(H::HMatrix) = (H.d, H.d)

	Base.getindex(A::HMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*(i==j) + one(T)
	Base.getindex(::HMatrix{T}, i::Integer, j::Integer) where {T} = one(T)*(i == j) + one(T)

Add elementary matrices constructors #21

Add elementary matrices constructors #21

Conversation

rlschuller commented Sep 12, 2021

juliohm 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

rlschuller Sep 12, 2021 • 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

rlschuller commented Sep 12, 2021

juliohm commented Sep 12, 2021 • edited Loading

codecov-commenter commented Sep 12, 2021 • edited Loading

Codecov Report

juliohm commented Sep 12, 2021 via email

juliohm commented Sep 12, 2021

rlschuller Sep 12, 2021 •

edited

Loading

juliohm commented Sep 12, 2021 •

edited

Loading

codecov-commenter commented Sep 12, 2021 •

edited

Loading