JuliaLang / julia

The Julia Programming Language
https://julialang.org/
MIT License
45.61k stars 5.48k forks source link

Constants Zeros and Ones that behave somewhat like I #46477

Open ojwoodford opened 2 years ago

ojwoodford commented 2 years ago

In scientific papers using matrices, you'll often see notation where block matrices consist of blocks such as A (a full matrix), but also I (an appropriately sized identity matrix), 0 and 1 (appropriately sized blocks of all zeros or all ones).

Julia (in particular the LinearAlgebra standard library) has a nice constant, I, that creates an implicit identity matrix, so we can do things like:

D = [A I; I B] * C

and Julia will magically do the right thing. This makes it really easy to transcribe matrix algebra using the identity matrix into code, and makes it easy for scientists reading the code to understand what it does. And of course it is performant because the code understands what the identity is, so doesn't waste time doing unnecessary multiplications by 0 and 1.

It would be really great to have additional constants, Zeros and Ones, which represent blocks of zeros and ones respectively, so we can similarly write:

D = [A Zeros; Ones B] * C

and again Julia will magically do the right thing, again not wasting time with unnecessary multiplications. This way transcribing matrix algebra from papers becomes even easier.

Of course, I is always assumed to be a square matrix, whereas 0 and 1 are not, so sometimes the dimensionality of the block cannot be inferred. In cases such as:

C = [A Zeros; Ones B]
C = A * [B; Ones] 

it can be inferred. However, in cases such as:

C = [A Zeros]
C = [A Ones] * [B; Ones] 

one of the dimensions of each block cannot be inferred. In this case, a compile-time or run-time block size parameter can be specified for the unknown dimension, e.g. Ones{1} or Ones(1). This could default to one, such that [A; Ones] would implicitly concatenate a row of ones to the bottom of A.

jishnub commented 2 years ago

If you're only working with square matrices, 0I should give you a block of zeros. In general, this sounds like something FillArrays may be able to implement, given that they have Zeros and Ones as array types already, albeit with the sizes specified explicitly. This use case would require one to figure the sizes out automatically.

mcabbott commented 2 years ago

Matrix(I, 3, 4) allows non-square I. It's a bit unfortunate that A = rand(3,4); [A I; I A] already does something, but rows not blocks, it's like vcat([A I], [I A]).