Re-posted from: https://blog.kdheepak.com/advent-of-code-2020-retrospective/index.html
It is Christmas night, and it is the first time this month that I haven’t had to plan my schedule for the evening around a programming puzzle contest. For the last 25 days this month, I participated in Advent of Code 2020, and I managed to collect all 50 stars!
Re-posted from: https://bkamins.github.io/julialang/2020/12/24/minilanguage.html
As it is end-of-year I thought of writing a longer post that would be a useful
reference to DataFrames.jl users.
In DataFrames.jl we have five functions that can be used to perform
transformations of columns of a data frame:
combine: create a new DataFrame populated with columns that are results ofselect: create a new DataFrame that has the same number of rows as theselect(df) which produces an empty data frame);select!: the same as select but updates the passed data frame in place;transform: the same as select but keeps the columns that were alreadytransform);transform!: the same as transform but updates the passed data frame inThe same functions also work with GroupedDataFrame with the difference
that the transformations are applied to groups and then combined. Here, an
important distinction is that combine again allows transformations to produce
any number of rows and they are combined in the order of groups in the
GroupedDataFrame. On the other hand select, select!, transform and
transform! require transformations to produce the same number of rows for each
group as the source group and produce a result that has the same row order
as the parent data frame of GroupedDataFrame passed. This rule has two
important implications:
select, select!, transform andtransform! operations on a GroupedDataFrame whose groups do not cover allparent data frame;select, select!, transform and transform!, contrary to combine,GroupedDataFrame.In this post I want to explain what I mean when I write transformations.
These transformations follow a so-called DataFrames.jl minilanguage that is
largely based on using => operator. You can find a specification of this
minilanguage here. However, because the rules have to be precise they
are relatively hard to read. Therefore in this post I will introduce them by
example. I will give all the examples with data frames as a source, but as noted
above they naturally extend to GroupedDataFrame case so I will add the
examples for GroupedDataFrame only in cases when it is especially relevant (in
order to read these parts of the post please earlier check out how groupby
function works in DataFrames.jl).
This post is written under Julia 1.5.3 and DataFrames.jl 0.22.2.
The data frame that we are going to use in our examples is defined follows
(I sampled some random data to populate it, with the exception of one name
to celebrate this recent masterpiece):
julia> df = DataFrame(id = 1:6,
name = ["Aaron Aardvark", "Belen Barboza",
"春 陈", "Даниил Дубов",
"Elżbieta Elbląg", "Felipe Fittipaldi"],
age = [50, 45, 40, 35, 30, 25],
eye = ["blue", "brown", "hazel", "blue", "green", "brown"],
grade_1 = [95, 90, 85, 90, 95, 90],
grade_2 = [75, 80, 65, 90, 75, 95],
grade_3 = [85, 85, 90, 85, 80, 85])
6×7 DataFrame
Row │ id name age eye grade_1 grade_2 grade_3
│ Int64 String Int64 String Int64 Int64 Int64
─────┼────────────────────────────────────────────────────────────────────
1 │ 1 Aaron Aardvark 50 blue 95 75 85
2 │ 2 Belen Barboza 45 brown 90 80 85
3 │ 3 春 陈 40 hazel 85 65 90
4 │ 4 Даниил Дубов 35 blue 90 90 85
5 │ 5 Elżbieta Elbląg 30 green 95 75 80
6 │ 6 Felipe Fittipaldi 25 brown 90 95 85
These are the most simple operations that the discussed functions can perform.
In this case the result of combine and select are the same (using
transform is not very useful here, as it retains all the columns anyway).
The allowed column selectors are: column number, column name, a regular
expression, Not, Between, Cols, All(), and :. In order to select a
column or a set of columns you just pass them as arguments. Here is an example
of a simple column selection:
julia> select(df, r"grade", :)
6×7 DataFrame
Row │ grade_1 grade_2 grade_3 id name age eye
│ Int64 Int64 Int64 Int64 String Int64 String
─────┼────────────────────────────────────────────────────────────────────
1 │ 95 75 85 1 Aaron Aardvark 50 blue
2 │ 90 80 85 2 Belen Barboza 45 brown
3 │ 85 65 90 3 春 陈 40 hazel
4 │ 90 90 85 4 Даниил Дубов 35 blue
5 │ 95 75 80 5 Elżbieta Elbląg 30 green
6 │ 90 95 85 6 Felipe Fittipaldi 25 brown
Above, we use a common pattern showing how one can move columns to the front of
the data frame. The r"grade" regular expression picks all columns that contain
"grade" string in their name, then : picks all the remaining columns.
An important rule here is that if we pass several column selectors that pick
multiple columns, as you can see above, it is allowed that they select the same
columns and they are included in the target data frame using first-in-first-out
policy. However, it is not allowed to select the same column twice using single
column selectors, so this works:
julia> select(df, "name", [2])
6×1 DataFrame
Row │ name
│ String
─────┼───────────────────
1 │ Aaron Aardvark
2 │ Belen Barboza
3 │ 春 陈
4 │ Даниил Дубов
5 │ Elżbieta Elbląg
6 │ Felipe Fittipaldi
but this fails:
julia> select(df, "name", 2)
ERROR: ArgumentError: duplicate output column name: :name
The same rule applies also to more advanced transformations that I cover below
(i.e. duplicates are allowed only in plain column selectors picking multiple
columns).
The simplest way to specify a transformation is
source_column => transformation => target_column_name
In this scenario the source_column is passed as an argument to
transformation and stored in target_column_name column.
Here is an example:
julia> using Statistics
julia> select(df, :age => mean => :meanage)
6×1 DataFrame
Row │ meanage
│ Float64
─────┼─────────
1 │ 37.5
2 │ 37.5
3 │ 37.5
4 │ 37.5
5 │ 37.5
6 │ 37.5
julia> combine(df, :age => mean => :meanage)
1×1 DataFrame
Row │ meanage
│ Float64
─────┼─────────
1 │ 37.5
Observe that because select produces as many rows in the produced data frame
as there are rows in the source data frame, a single value is repeated
accordingly. This is not the case for combine. However, if other columns
in combine would produce multiple rows the repetition also happens:
julia> combine(df, :age => mean => :meanage, :eye => unique => :eye)
4×2 DataFrame
Row │ meanage eye
│ Float64 String
─────┼─────────────────
1 │ 37.5 blue
2 │ 37.5 brown
3 │ 37.5 hazel
4 │ 37.5 green
Note, however, that it is not allowed to return vectors of different lengths in
different transformations:
julia> combine(df, :age, :eye => unique => :eye)
ERROR: ArgumentError: New columns must have the same length as old columns
Let us discuss one more example, this case using GroupedDataFrame. As you
can see below vectors get expanded into multiple columns by default, e.g.:
julia> combine(groupby(df, :eye),
:name => (x -> identity(x)) => :name,
:age => mean => :age)
6×3 DataFrame
Row │ eye name age
│ String String Float64
─────┼────────────────────────────────────
1 │ blue Aaron Aardvark 42.5
2 │ blue Даниил Дубов 42.5
3 │ brown Belen Barboza 35.0
4 │ brown Felipe Fittipaldi 35.0
5 │ hazel 春 陈 40.0
6 │ green Elżbieta Elbląg 30.0
(as a side note observe that I used (x -> identity(x)) transformation although
in this case it would be enough just to write :name; the reason is that I
wanted to highlight that if you use an anonymous function in the transformation
definition you have to wrap it in ( and ) as shown above; otherwise the
expression is parsed by Julia in an unexpected way due to the => operator
precedence level)
However, in some cases it would be more natural not to expand the :name
column into multiple rows. You can easily achieve it by protecting the value
returned by the transformation using Ref (just as in broadcasting):
julia> combine(groupby(df, :eye),
:name => (x -> Ref(identity(x))) => :name,
:age => mean => :age)
4×3 DataFrame
Row │ eye name age
│ String SubArray… Float64
─────┼────────────────────────────────────────────────────
1 │ blue ["Aaron Aardvark", "Даниил Дубов… 42.5
2 │ brown ["Belen Barboza", "Felipe Fittip… 35.0
3 │ hazel ["春 陈"] 40.0
4 │ green ["Elżbieta Elbląg"] 30.0
In this case the vectors produced by our transformation are kept in a single row
of the produced data frame. Note that in this case we also could have just
written:
julia> combine(groupby(df, :eye), :name => Ref => :name, :age => mean => :age)
4×3 DataFrame
Row │ eye name age
│ String SubArray… Float64
─────┼────────────────────────────────────────────────────
1 │ blue ["Aaron Aardvark", "Даниил Дубов… 42.5
2 │ brown ["Belen Barboza", "Felipe Fittip… 35.0
3 │ hazel ["春 陈"] 40.0
4 │ green ["Elżbieta Elbląg"] 30.0
as we were applying identity transformation to our data, which can be skipped.
Often you want to apply some function not to the whole column of a data frame,
but rather to its individual elements. Normally one would achieve this using
broadcasting like this:
julia> select(df, :name => (x -> uppercase.(x)) => :NAME)
6×1 DataFrame
Row │ NAME
│ String
─────┼───────────────────
1 │ AARON AARDVARK
2 │ BELEN BARBOZA
3 │ 春 陈
4 │ ДАНИИЛ ДУБОВ
5 │ ELŻBIETA ELBLĄG
6 │ FELIPE FITTIPALDI
This pattern is encountered very often in practice, therefore there is a ByRow
convenience wrapper for a function that creates its broadcasted variant:
julia> select(df, :name => ByRow(uppercase) => :NAME)
6×1 DataFrame
Row │ NAME
│ String
─────┼───────────────────
1 │ AARON AARDVARK
2 │ BELEN BARBOZA
3 │ 春 陈
4 │ ДАНИИЛ ДУБОВ
5 │ ELŻBIETA ELBLĄG
6 │ FELIPE FITTIPALDI
One can skip specifying a target column name, in which case it is generated
automatically by suffixing source column name by function name that is applied
to it, e.g.:
julia> select(df, :name => ByRow(uppercase))
6×1 DataFrame
Row │ name_uppercase
│ String
─────┼───────────────────
1 │ AARON AARDVARK
2 │ BELEN BARBOZA
3 │ 春 陈
4 │ ДАНИИЛ ДУБОВ
5 │ ELŻBIETA ELBLĄG
6 │ FELIPE FITTIPALDI
If you want to avoid this then pass renemecols=false keyword argument (this
is mostly useful if you perform transformations):
julia> select(df,
:grade_1 => ByRow(x -> x / 100),
:grade_2 => ByRow(x -> x / 100),
:grade_3 => ByRow(x -> x / 100),
renamecols=false)
6×3 DataFrame
Row │ grade_1 grade_2 grade_3
│ Float64 Float64 Float64
─────┼───────────────────────────
1 │ 0.95 0.75 0.85
2 │ 0.9 0.8 0.85
3 │ 0.85 0.65 0.9
4 │ 0.9 0.9 0.85
5 │ 0.95 0.75 0.8
6 │ 0.9 0.95 0.85
Similarly to skipping target column name you can skip the transformation part of
the pattern we discussed, in which case it is assumed to be the identity
function. In effect we get a way to rename columns in transformations:
julia> select(df, r"grade", :grade_1 => :g1, :grade_2 => :g2, :grade_3 => :g3)
6×6 DataFrame
Row │ grade_1 grade_2 grade_3 g1 g2 g3
│ Int64 Int64 Int64 Int64 Int64 Int64
─────┼────────────────────────────────────────────────
1 │ 95 75 85 95 75 85
2 │ 90 80 85 90 80 85
3 │ 85 65 90 85 65 90
4 │ 90 90 85 90 90 85
5 │ 95 75 80 95 75 80
6 │ 90 95 85 90 95 85
If you want to perform multiple column transformations that have a similar
structure (as in the example above) you can pass a vector of transformations:
julia> select(df,
[x => ByRow(x -> x / 100) for x in [:grade_1, :grade_2, :grade_3]],
renamecols=false)
6×3 DataFrame
Row │ grade_1 grade_2 grade_3
│ Float64 Float64 Float64
─────┼───────────────────────────
1 │ 0.95 0.75 0.85
2 │ 0.9 0.8 0.85
3 │ 0.85 0.65 0.9
4 │ 0.9 0.9 0.85
5 │ 0.95 0.75 0.8
6 │ 0.9 0.95 0.85
or even shorter using broadcasting and taking advantage of names function:
julia> select(df, names(df, r"grade") .=> ByRow(x -> x / 100), renamecols=false)
6×3 DataFrame
Row │ grade_1 grade_2 grade_3
│ Float64 Float64 Float64
─────┼───────────────────────────
1 │ 0.95 0.75 0.85
2 │ 0.9 0.8 0.85
3 │ 0.85 0.65 0.9
4 │ 0.9 0.9 0.85
5 │ 0.95 0.75 0.8
6 │ 0.9 0.95 0.85
You now know almost all about single column selection. The only thing left to
learn is that nrow function has a special treatment. It does not require
passing source column name (but allows specifying of target column name; the
default name is :nrow) and produces number of rows in the passed data frame.
Here is an example of both:
julia> combine(groupby(df, :eye), nrow, nrow => :count)
4×3 DataFrame
Row │ eye nrow count
│ String Int64 Int64
─────┼──────────────────────
1 │ blue 2 2
2 │ brown 2 2
3 │ hazel 1 1
4 │ green 1 1
That was a lot of information. Probably this is a good moment to take a short
break.
Next we move on to the cases when there are multiple source columns or multiple
columns are produced as a result of the transformation.
The simplest way to pass multiple columns to a transformation is to pass their
list in a vector. In this case these columns are passed as consecutive
positional arguments to the function. Here is an example (also observe how
automatic column naming is done in this case):
julia> combine(df,
[:age, :grade_1] => cor,
[:age, :grade_2] => cor,
[:age, :grade_3] => cor)
1×3 DataFrame
Row │ age_grade_1_cor age_grade_2_cor age_grade_3_cor
│ Float64 Float64 Float64
─────┼───────────────────────────────────────────────────
1 │ 0.0710072 -0.536745 0.338062
Alternatively you can pass the columns as a single argument being a
NamedTuple, in order to achieve this wrap the list of passed columns in
AsTable. Here is a simple example:
julia> combine(df, :grade_1, :grade_2,
AsTable([:grade_1, :grade_2]) => ByRow(x -> x.grade_1 > x.grade_2))
6×3 DataFrame
Row │ grade_1 grade_2 grade_1_grade_2_function
│ Int64 Int64 Bool
─────┼────────────────────────────────────────────
1 │ 95 75 true
2 │ 90 80 true
3 │ 85 65 true
4 │ 90 90 false
5 │ 95 75 true
6 │ 90 95 false
Let us see these two rules at play in a common task of finding a sum of several
columns:
julia> select(df, names(df, r"grade") => +, AsTable(names(df, r"grade")) => sum)
6×2 DataFrame
Row │ grade_1_grade_2_grade_3_+ grade_1_grade_2_grade_3_sum
│ Int64 Int64
─────┼────────────────────────────────────────────────────────
1 │ 255 255
2 │ 255 255
3 │ 240 240
4 │ 265 265
5 │ 250 250
6 │ 270 270
Normally all the transformation functions assume that a single column is
returned from a transformation function. Here is a more advanced example:
julia> select(df,
:name => ByRow(x -> (; ([:firsname, :lastname] .=> split(x))...)))
6×1 DataFrame
Row │ name_function
│ NamedTuple…
─────┼───────────────────────────────────
1 │ (firsname = "Aaron", lastname = …
2 │ (firsname = "Belen", lastname = …
3 │ (firsname = "春", lastname = "陈…
4 │ (firsname = "Даниил", lastname =…
5 │ (firsname = "Elżbieta", lastname…
6 │ (firsname = "Felipe", lastname =…
In the above code we have used the following pattern that is a convenient way
of programmatic creation of NamedTuples:
julia> (; :a => 1, :b => 2)
(a = 1, b = 2)
julia> (; [:a => 1, :b => 2]...)
(a = 1, b = 2)
However, it is natural to ask if it is possible to produce the separate
:firstname and :lastname columns in the data frame. You can do it by
adding AsTable as target column names specifier:
julia> select(df, :name =>
ByRow(x -> (; ([:firsname, :lastname] .=> split(x))...)) =>
AsTable)
6×2 DataFrame
Row │ firsname lastname
│ SubStrin… SubStrin…
─────┼───────────────────────
1 │ Aaron Aardvark
2 │ Belen Barboza
3 │ 春 陈
4 │ Даниил Дубов
5 │ Elżbieta Elbląg
6 │ Felipe Fittipaldi
In general it is not required that one has to produce a vector of NamedTuples
to get this result. It is enough to have any standard iterable (e.g. a vector
or a Tuple):
julia> select(df, :name => ByRow(split) => AsTable)
6×2 DataFrame
Row │ x1 x2
│ SubStrin… SubStrin…
─────┼───────────────────────
1 │ Aaron Aardvark
2 │ Belen Barboza
3 │ 春 陈
4 │ Даниил Дубов
5 │ Elżbieta Elbląg
6 │ Felipe Fittipaldi
Note that as split produces a vector (without column names) the names are
generated automatically as :x1 and :x2. You can override this behavior by
specifying the target column names explicitly:
julia> select(df, :name => ByRow(split) => [:firsname, :lastname])
6×2 DataFrame
Row │ firsname lastname
│ SubStrin… SubStrin…
─────┼───────────────────────
1 │ Aaron Aardvark
2 │ Belen Barboza
3 │ 春 陈
4 │ Даниил Дубов
5 │ Elżbieta Elbląg
6 │ Felipe Fittipaldi
Several values that can be returned by a transformation are treated to produce
multiple columns by default. Therefore they are not allowed to be returned
from a function unless AsTable or multiple target column names are specified.
Here is an example:
julia> combine(df, :age => x -> (age=x, age2 = x.^2))
ERROR: ArgumentError: Table returned but a single output column was expected
julia> combine(df, :age => (x -> (age=x, age2 = x.^2)) => AsTable)
6×2 DataFrame
Row │ age age2
│ Int64 Int64
─────┼──────────────
1 │ 50 2500
2 │ 45 2025
3 │ 40 1600
4 │ 35 1225
5 │ 30 900
6 │ 25 625
The list of return value types treated in this way consists of four elements:
AbstractDataFrame,DataFrameRow,AbstractMatrix,NamedTuple.In the example above we returned a NamedTuple. Note, however, that returning
a NamedTuple is not the same as returning a vector of NamedTuples (as we
did in the topmost example in this section) as this is allowed.
For the purposes of detection of column names conflicts multiple columns
returned are treated as a series of single returned columns (so duplicates are
not allowed):
julia> select(df, :name => ByRow(split) => [:name, :lastname], :name)
ERROR: ArgumentError: duplicate output column name: :name
Just for completeness of discussion note that this works:
julia> select(df, :name => ByRow(split) => [:name, :lastname], [:name])
6×2 DataFrame
Row │ name lastname
│ SubStrin… SubStrin…
─────┼───────────────────────
1 │ Aaron Aardvark
2 │ Belen Barboza
3 │ 春 陈
4 │ Даниил Дубов
5 │ Elżbieta Elbląg
6 │ Felipe Fittipaldi
as [:names] is considered to be a multiple column selector (thus it is not
producing a duplicate and is just ignored in this case).
The traditional transformation style in the combine function was to pass a
transformation as a function (without using the => minilanguage). This style
is also allowed currently in all transformation functions. Additionally if you
want to use this approach you can pass the function as a first argument of the
function (which is often convenient in combination with do block style).
In this traditional style we are not allowed to specify source columns.
Therefore such a function always takes a data frame (it is most useful with
GroupedDataFrame where the data frame is representing a given group of the
source data frame)
Similarly traditional style does not allow specifying target column names.
Therefore the following defaults are taken:
AbstractDataFrame, DataFrameRow, AbstractMatrix, or NamedTuple isNamedTuple it must contain only vectors or only single values, as if theyHere are examples of this style:
julia> combine(groupby(df, :eye)) do sdf
return mean(sdf.age)
end
4×2 DataFrame
Row │ eye x1
│ String Float64
─────┼─────────────────
1 │ blue 42.5
2 │ brown 35.0
3 │ hazel 40.0
4 │ green 30.0
julia> combine(groupby(df, :eye), sdf -> mean(sdf.age))
4×2 DataFrame
Row │ eye x1
│ String Float64
─────┼─────────────────
1 │ blue 42.5
2 │ brown 35.0
3 │ hazel 40.0
4 │ green 30.0
julia> combine(groupby(df, :eye),
sdf -> (count = nrow(sdf),
mean_age = mean(sdf.age),
std_age = std(sdf.age)))
4×4 DataFrame
Row │ eye count mean_age std_age
│ String Int64 Float64 Float64
─────┼───────────────────────────────────
1 │ blue 2 42.5 10.6066
2 │ brown 2 35.0 14.1421
3 │ hazel 1 40.0 NaN
4 │ green 1 30.0 NaN
As a special rule for processing of GroupedDataFrame object is that
it is not allow to change the values stored in the grouping columns as they
are kept by default.
Here is a simple example:
julia> combine(groupby(df, :eye), :eye => ByRow(uppercase) => :eye)
ERROR: ArgumentError: column :eye in returned data frame is not equal to grouping key :eye
However, you are allowed to do this if you drop the grouping columns by passing
keepkeys=false keyword argument:
julia> transform(groupby(df, :eye),
:eye => ByRow(uppercase) => :eye,
keepkeys=false)
6×7 DataFrame
Row │ id name age eye grade_1 grade_2 grade_3
│ Int64 String Int64 String Int64 Int64 Int64
─────┼────────────────────────────────────────────────────────────────────
1 │ 1 Aaron Aardvark 50 BLUE 95 75 85
2 │ 2 Belen Barboza 45 BROWN 90 80 85
3 │ 3 春 陈 40 HAZEL 85 65 90
4 │ 4 Даниил Дубов 35 BLUE 90 90 85
5 │ 5 Elżbieta Elbląg 30 GREEN 95 75 80
6 │ 6 Felipe Fittipaldi 25 BROWN 90 95 85
(note that for transform the key columns would be retained in the produced
data frame even with keepkeys=false; in this case this keyword argument only
influences the fact if we check that key columns have not changed in this case)
This was long. As a conclusion let me comment that all the above-mentioned
styles can be freely combined in a single transformation call:
julia> combine(groupby(df, :eye),
r"grade", names(df, r"grade") => (+) => :total_grade,
sdf -> nrow(sdf), nrow, :age => :AGE)
6×8 DataFrame
Row │ eye grade_1 grade_2 grade_3 total_grade x1 nrow AGE
│ String Int64 Int64 Int64 Int64 Int64 Int64 Int64
─────┼─────────────────────────────────────────────────────────────────────
1 │ blue 95 75 85 255 2 2 50
2 │ blue 90 90 85 265 2 2 35
3 │ brown 90 80 85 255 2 2 45
4 │ brown 90 95 85 270 2 2 25
5 │ hazel 85 65 90 240 1 1 40
6 │ green 95 75 80 250 1 1 30
If you got here then congratulations – you have mastered the DataFrames.jl
transformation minilanguage.
By: Julia on μβ
Re-posted from: https://matbesancon.github.io/post/2020-12-24-chains_sets2/
In a previous post, I detailed some of the features of
MathOptSetDistances.jl
and the evolution of the idea behind it. This is part II focusing on derivatives.
The most interesting part of the packages is the projection onto a set.
For some applications, what we need is not only the projection but also the
derivative of this projection.
One answer here would be to let Automatic Differentiation (AD) do the work.
However:
For these reasons, AD systems often let users implement some derivatives themselves,
but as a library developer, I do not want to depend on a full AD package
(and force downstream users to do so).
ChainRules.jl is a Julia package
addressing exactly the issue mentioned above: it defines a set of primitives
to talk about derivatives in Julia.
Library developers can implement custom derivatives for their own functions and types.
Finally, AD library developers can leverage ChainRules.jl to obtain derivatives
from functions when available, and otherwise use AD mechanisms to obtain them from
more elementary functions.
The logic and motivation is explained in more details in Lyndon’s talk
at JuliaCon 2020 and the package documentation
which is very instructive on AD in general.
We are interested in computing
$D\Pi_{\mathcal{S}}(v)$, the derivative of the projection with respect to the
initial point. As a refresher, if $\Pi_s(\cdot)$ is a function from $V$ onto itself,
and if $V$ then the derivative $D\Pi$ maps a point in $V$ onto a linear map
from the *tangent space* of $V$ onto itself.
The tangent space of $V$ is roughly speaking the space where differences of
values in $V$ live. If $V$ corresponds to real numbers, then the tangent space
will also be real numbers, but if $V$ is a space of time/dates, then the tangent
space is a duration/time period. See here3 for more references.
Again, roughly speaking, this linear map takes perturbations of the input $\Delta v$
and maps them to perturbation of the projected point $\Delta v_p$.
As an example warm-up:
$D\Pi_{\mathcal{S}}(v)$ is a linear map from $\mathcal{V}$ to $\mathcal{V}$.
If $v \in \mathbb{R}^n$, it can be represented as a
$n\times n$ matrix.
There are several ways of representing linear maps, see the LinearOperators.jl
package for some insight. Two approaches (for now) are implemented for set distances:
(1) has been implemented by Akshay for many sets
during his GSoC this summer, along with the projections themselves.
(1) corresponds to computing the derivative eagerly as a full matrix, thus
paying storage and computation cost upfront. The advantage is the simplicity for standard vectors,
take v, s, build and return the matrix.
(2) is the building block for forward-mode differentiation:
given a point $v$ and an input perturbation $\Delta v$, compute the output perturbation.
(3) corresponds to a building block for reverse-mode differentiation.
An aspect of the matrix approach is that it works well for 1-D arrays
but gets complex quite quickly for other structures, including multi-argument
functions or matrices. Concatenating everything into a vector is too rigid.
The nonnegative orthant cone is the set $\mathbb{R}^n_+$; it is represented in MOI
as MOI.Nonnegatives(n) with n the dimension.
The projection is simple because it can be done elementwise:
$$
(\Pi_S(v))_i = max(v_i, 0) \,\,\forall i.
$$
In other terms, any non-diagonal term of the gradient matrix is 0 for any $v$.
Here is a visualization made with haste for $n=2$ using the very promising Javis.jl:
The red circle is a vector in the plane and the blue square its projection.4
The Julia implementation follows the same idea, here in a simplified version:
function projection_on_set(v::AbstractVector{T}, s::MOI.Nonnegatives) where {T}
return max.(v, zero(T))
end
For each component $i \in 1..n$, there are two cases to compute its derivative, either
the constraint is active or not.
$$
\begin{align}
v_i < 0 & \Rightarrow \frac{\partial \Pi_i}{\partial v_i}(v) = 0\\
v_i > 0 & \Rightarrow \frac{\partial \Pi_i}{\partial v_i}(v) = 1.
\end{align}
$$
The projection is not differentiable on points where one of the components is 0.
The convention usually taken is to return any quantity on such point
(to the best of my knowledge, no system guarantees a subgradient).
The Julia implementation holds on two lines:
function projection_gradient_on_set(v::AbstractVector{T}, ::MOI.Nonnegatives) where {T}
y = (sign.(v) .+ one(T)) / 2
return LinearAlgebra.Diagonal(y)
end
First the diagonal of the matrix is computed using broadcasting and the sign function.
Then a LinearAlgebra.Diagonal matrix is constructed. This matrix type is sparsity-aware,
in the sense that it encodes the information of having only non-zero entries on
the diagonal. We save on space, using $O(n)$ memory instead of $O(n^2)$ for a
full matrix, and can benefit from specialized methods down the line.
We implemented the matrix approach from scratch. Even though we materialize the
derivative as a diagonal matrix, it still costs storage, which will become a
burden when we compose this projection with other functions and compute derivatives
on the composition.
For a function f, value v and tangent Δv, the forward rule, or frule
in ChainRules.jl does two things at once:
y = f(v),∂y = Df(v) Δv.The motivation for computing the two values at once is detailed in the
documentation.
Quite often, computing the derivative will require computing f(v) itself
so it is likely to be interesting to return it anyway instead of forcing the user
to call the function again.
The exact signature of ChainRulesCore.frule involves some details we want to
ignore for now, but the essence is as follows:
function frule((Δself, v...), ::typeof(f), v...; kwargs...)
...
return y, ∂y
end
∂Y is the directional derivative using the direction Δx. Note here the variadic
Δx and x, since we do not want to impose a rigid, single-argument structure
to functions. The Δself argument is out of scope for this post but you can read
on its use in the docs.
For our set projection, it may look like this:
function ChainRulesCore.frule(
(_, Δv, _),
::typeof(projection_on_set),
v::AbstractVector{T}, s::MOI.Nonnegatives) where {T}
vproj = projection_on_set(v, s)
∂vproj = Δv .* (v .>= 0)
return vproj, ∂vproj
end
The last computation line leverages broadcast to express elementwise the
multiplication of Δv with the indicator of v[i] being nonnegative.
The important thing to note here is that we never build the derivative as a data
structure. Instead, we implement it as a function. An equivalent using our
projection_gradient_on_set would be:
function projection_directional_derivative(v, Δv, s)
vproj = projection_on_set(v, s)
DΠ = projection_gradient_on_set(v, s)
∂vproj = DΠ * Δv
return vproj, ∂vproj
end
Notice the additional allocation and matrix-vector product.
The forward mode is fairly intuitive, the backward mode less so.
The motivation for using it, and the reason it is the favoured one for several
important fields using AD, is that it can differentiate a composition of functions
with only matrix-vector products, instead of requiring matrix-matrix products.
What it computed is, given a perturbation in the output (or seed), provide the
corresponding perturbation in the input.
There are great resources online which will explain it in better terms than I could
so we will leave it at that.
Looking at the rrule signature from ChainRules.jl:
function rrule(::typeof(f), x...; kwargs...)
y = f(x...)
function pullback_f(Δy)
# implement the pullback here
return ∂self, ∂x
end
return y, pullback_f
end
This is a bit denser. rrule takes the function as input and its arguments.
So far so good. It returns two things, the value y of the function, similalry to frule
and a pullback. This term comes from differential geometry and in the context
of AD, is also referred to as a backpropagator. Again, the ChainRules
docs
got your back with great explanations.
It also corresponds to the Jacobian-transpose vector product if you prefer the term.
In the body of pullback_f, we compute the variation of the output with respect to each input.
If we give the pullback a 1 or 1-like as input, we compute the gradient,
the partial derivative of f with respect to each input x[i] evaluated at the
point x.
Here is the result for our positive orthant (again, simplified for conciseness):
function ChainRulesCore.rrule(::typeof(projection_on_set), v, s::MOI.Nonnegatives)
vproj = projection_on_set(v, s)
function pullback(Δvproj)
n = length(v)
v̄ = zeros(eltype(Δvproj), n)
for i in 1:n
if vproj[i] == v[i]
v̄[i] = Δvproj[i]
end
end
return (ChainRulesCore.NO_FIELDS, v̄, ChainRulesCore.DoesNotExist())
end
return (vproj, pullback)
end
The first step is computing the projection, here we do not bother with saving
for loops and just call the projection function.
For each index i of the vector, if the i-th projection component is equal to
the i-th initial point, $v_i$ is in the positive orthant and variations of
the output are directly equal to variations of the input. Otherwise,
this means the non-negativity constraint is tight, the projection lies on
the boundary vproj[i] = 0, and output variations are not propagated to the input
since the partial derivative is zero.
We see here that a tuple of 3 elements is returned. The first corresponds to
∂self, out of the scope for this package. The second is the interesting one,
v̄, the derivative with respect to the input point.
The last one ChainRulesCore.DoesNotExist() indicates that there is no derivative
with respect to the last argument of projection_on_set, namely the set s.
This makes sense because there is nothing to differentiate in the set.
An interesting point to notice is that the implementation, not the types defines the derivatives.
A non-trivial example would be a floating-point argument p only used to extract
the sign bit. This means it would not have a notion of local perturbation.
The type (a floating-point) would be interpreted as differentiable.
To my understanding, Swift for Tensorflow uses
a type-first approach, where types indicate what field gets differentiated.
If you imagine using this in practice, in an AD library for instance,
one would first call rrule forward, computing primal values and collecting the
successive pullbacks. Once we arrive at the end of our chain of functions,
we could backpropagate from $\Delta Y_{final} = 1$, walking our way back to
the primary input parameters.
This post comes after a few weeks of work on MathOptSetDistances.jl,
the package with the actual implementation of the presented features.
There is still a lot to learn and do on the topic, including solutions to more
projections and derivatives thereof, but also interesting things to build upon.
Defining derivatives and projections is after all a foundation for greater things to
happen.
See H. Friberg’s talk on exponential cone projection in Mosek at ISMP 2018 ↩︎
An example case for the projection onto the Positive Semidefinite cone ↩︎
If like me you haven’t spent much time lying around differential geometry books,
the ChainRules.jl
documentation has a great developer-oriented explanation.
For more visual explanations, Keno Fischer had a recent talk on
the topic. ↩︎