mcabbott

Covariant and contravariant indices would be the formal terms. I'm not really sure whether I've seen "upstairs" written down.

Sub/superscript... strike me as the typographical terms, not the meaning? Like $x_\mathrm{alice}$ is certainly a subscript, and footnote 2 is a superscript, but neither is an index.

This doesn't really help with programming, but in physics it's traditional to use up- and down-stairs indices, which makes the distinction you want very clear.

If input x has components xⁿ, and output f(x) components fᵐ, then the Jacobian is ∂ₙfᵐ which has one index upstairs and one downstairs. The derivative has a downstairs index... because x is in the denominator of d/dx, roughly? If x had units seconds, then d/dx has units per second.

Whereas if g(x) is a number, the gradient is ∂ₙg, and the Hessian is ∂ₙ∂ₙ₂g with two downstairs indices. You might call this a (0,2) tensor, while the Jaconian is (1,1). Most of the matrices in ordinary linear algebra are (1,1) tensors.

That's what Reactant.jl is aiming to do: Take the LLVM code from Julia and pipe to to XLA, where it can benefit from all the investment which make Jax fast.

Same author as the more mature Enzyme.jl, which is AD done at the LLVM level.

Julia's flexibility makes it easy to write a simple AD system, but perhaps we're learning that it's the wrong level to write the program transformations needed for really efficient AD. Perhaps the same is true of "run this on a GPU, or 10" transformations.

This is also how it works in Julia, where macros digest notation for einsum-like operations before compile-time. In fact the linked file's explanatory comment:

     (einsum (A i j) (B i k) (C k j)) 
    results in the the updates
      A[i,j] = \\sum_k B[i,k]C[k,j],
    which is equivalent to matrix multiplication.

very nearly contains the syntax used by all the Julia packages (where @ marks a macro), which is

    @tensor A[i,j] = B[i,k] * C[k,j]

(using https://github.com/Jutho/TensorOperations.jl, but see also OMEinsum, Finch, and my Tullio, TensorCast.)

    using Unitful: m; foo(2.0m)

With the above definition, this will give DimensionError: 2.0 m and 1.0 are not dimensionally compatible.

Probably this means you should define `foo(x::Number) = x + oneunit(x)` to respect the Number interface. But this interface isn't very strictly defined. I believe `Base.oneunit` was added after someone started writing the Unitful package -- building something useful in a legal grey area, and formalising later?

Yes, it obeys enough algebraic laws that calling it a derivative is useful.

But I don't think there's any underlying notion of small changes to something continuous. It is not the slope of some smooth function.

Or perhaps better "foreigner", legal alien, someone from over the border.

The alien derivative is an operation which moves you from one sector to another. A bit like crossing a branch cut.

I presume so.

Note that this pseudo-code is very nearly valid Julia code (with my package):

    @tullio (min) C[i,j] := A[i,k] + B[k,j]

Yes, although that seems like the easy half of this, making sure `struct NewNum <: AbstractFloat` defines everything. There aren't yet tools for this but they are easy to imagine. And missing methods do give errors.

The hard half seems to be correctness of functions which accept quite generic objects. For example writing `f(x::Number)` in order to allow units, means you also allow quaternions, but many functions doing that will incorrectly assume numbers commute. (And not caring is, for 99% of these, the intention. But it's not encoded anywhere.) Less obviously, we can differentiate many things by passing dual numbers through `f(x::Real)`, but this tends to find edge cases nobody thought of. Right now if your algorithm branches on `if det(X) == 0` (or say a check that X is upper triangular) then it will sometimes give wrong answers. This one should be fixed soon, but I am sure there are other subtleties.

The unicode epsilon isn't in the public API, it's describing the 3rd positional argument.

This was added recently, and for some reason the PR (1840) didn't fix the docs, which is bad. The Optimisers.jl version has an explanation: https://fluxml.ai/Optimisers.jl/dev/api/#Optimisers.RMSProp

The struct's field name is `eta`, but this is an internal detail. Its constructor takes a positional-only argument, no public name.

The greek letter is used in the documentation. And the reason is that every optimiser's documentation links to the original paper, and tries to follow that. If the Adam paper calls the two decay rates β1, β2, then staying close to that seems the least confusing option.

> > Learning rate (η): Amount by which gradients are discounted before updating the weights.

> so this is already explicit to anyone who reads the documentation. The quibble in the post is about the named parameter.

As far as I can tell it's a documentation complaint. He has to remember "η" from the line with the signature, past the line "Gradient descent optimizer with learning rate η ...", and a heading "Parameters" until the line quoted which explains this in full.

He says this is the API, but that's inaccurate. The API being explained is that the first positional argument is the learning rate. It's not a keyword argument, so you cannot supply it by name. What variable names are used in the code is private, and in fact the struct's field name is `eta` so that you can access it without typing greek.

If this makes the top 10 list (even the top 10 list of documentation complaints) then Flux is doing OK. Especially the top 10 list of a guy with a PhD in a mathematical field. (From the sort of university which used to require students to know latin & greek, too.)

On the GPU, operations which look like matrix multiplication are indeed quite slow. As you say there are many tricks, and Tullio doesn't (right now) know them. For operations which are more like broadcasting it does much better.

On the CPU, the situation is much better. With some help from LoopVectorization.jl (which optimises micro-kernels) it will often beat OpenBLAS at matrix multiplication. The best-case scenario is an operation which would otherwise be permutedims plus matrix multiplication, for which it will often be several times faster, by fusing these.

The notation above is shared by some other packages. TensorOperations.jl is always decomposes to known kernels (including on the GPU) and OMEinsum.jl usually does so (with a fallback to loops), both more like einsum. TensorCast.jl is more like einops, just notation for writing reshape/permute/slice operations.

There is a discussion here:

https://discourse.julialang.org/t/rust-julia-comparison-post...

So far, a bug fix improves it 4700 min -> 150 min, changing the type stored -> 20 min, further optimisations -> 15 min.

Thanks for the references. We're trying to write up a paper about all this work, and (to us) it seems obvious that differential geometry is the right framework.

What reverse mode AD talks informally about as "gradients" are cotangents, and if you are constrained to a submanifold of another manifold (such as the space of all symmetric matrices, embedded in the space of all NxN matrices) then your cotangent is a projection of that in the larger space. This reduces to the obvious thing for symmetric matrices, as you say, but there are less obvious AD-relevant cases (like the space of orthogonal matrices, or the space of "ranges", vectors with uniformly spaced elements) where getting it right with your bare hands is tricky.

I don't see any of the relevant terminology in these links, but they do seem to be thinking about related problems. Perhaps they have re-invented the wheel? Amused to read that standard cookbooks seem to have faithfully reproduced the recipe for a square wheel, for decades.

Some points have simple answers:

* Support for using @copied s.x = 1 where s is an immutable struct.

I believe the standard package for this is this:

https://github.com/jw3126/Setfield.jl#usage

* Currently macros have to be valid Julia code,

There are two kinds of macros, and it sounds like you want `macro foo_str(x)` which accepts anything. But the tradeoff is that you have to parse it, rather than getting an expression tree.

https://docs.julialang.org/en/v1/manual/metaprogramming/#met...

The big Alcoa press is 50_000 tons, not lbs. There are slightly larger ones now, 60, 75ktons. (22 tons is, of course, force that any rail yard could produce a century ago.)

The Julia translations of these might look as follows:

    using OffsetArrays
    
    cash_flows = zeros(0:max_time)
    
    altitude = zeros(Bool, -100:100)
    
    for I in eachindex(cash_flows)
        some_function(cash_flows[I])
    end
    # or in this case avoiding indexing:
    for cf in cash_flows
        some_function(cf)
    end
    foreach(some_function, cash_flows)

OffsetArrays is not quite a standard library but it's close. A lot of library code will work like this, calling `eachindex` or `axes` so as to be indifferent to how the arrays are indexed, and to pass this behaviour through to outputs as appropriate.

The goal of Tullio.jl is certainly not to beat MKL at its own game. But it can still be 4-5 times faster on some permutedims + matmul operations, by fusing them.

Besides handling operations which aren't standard kernels, handling weird number types efficiently would be nice. (Being able to compile lighter-weight Julia images without BLAS libraries might be nice too, for other purposes.)

Sure! The desire for fancier notation is from more complicated examples, the kind where you write elaborate comments to explain why you're permuting 3rd & 4th dim of B to line up with C and part of D.

> until next time you have to read that

But how will this tool help? If you put the sample you used in a comment, then your confusion is wondering if that's in sync... and whether the input is like what you wrote.

IMO the sweet spot that's easy both for humans & computers is index notation. Instead of `tf.add(X, tf.expand_dims(Y, 1))` (their example, which I think also applies to 4-tensors in ways not trivial to see from the sample) something like "out_{r,c} = X_r + Y_c" is completely unambiguous, and quicker to type than making a sample.

I don't think you can do that in tf.einsum, nor einops. Is there a Python package which does this? In Julia you can write `@cast out[r,c] := X[r] + Y[c]` (with my package, but several others share this notation).

Right, NamedDims.jl is the one closest what the article describes. Very simple and light-weight, plays well with others.

There is a small zoo of packages attaching also labels along the indices, alla python's xarray. Perhaps it's a little too easy to write such things... but once we understand the design space the hope is to quietly take most of them to the woodshed.

You can do this:

    @_ images |> filter(size(_,1)>100, __) .|>
    resize(__, 100,100) .|>
    crop(__, 25,25,75,75)

with Underscores.jl. Not exactly `images.first.resize(100, 100).crop(25, 25, 75, 75)` anymore, to point out that this doesn't need every function to prioritise one argument, the pipe feeds into __ in any position. And that .|> applies them to each element.

At 2:39 it's saying that "doing good" and "doing well" are both correct, but mean different things: superman is doing good deeds, and scrooge is doing very well. But what to make of "I am good" isn't addressed (unless I missed it).

Those are all good synonyms for "got" in the context of shooting at things. But none of the others already has a strong meaning in the context of aircraft, and this other meaning does create some confusion, which is why many speakers would avoid it (if thinking clearly).

There are a bunch of packages exploring ideas around labelled dimensions, although we are still far from the vision of caring as little as we do about registers:

* https://github.com/JuliaArrays/AxisArrays.jl * https://github.com/davidavdav/NamedArrays.jl * https://github.com/invenia/NamedDims.jl * https://github.com/rafaqz/DimensionalData.jl * https://github.com/ITensor/ITensors.jl

While these don't pass dimension labels along in the type, there are also a number of Einsum-like packages which should work on any storage order, and could be taught to interface with such types:

* https://github.com/Jutho/TensorOperations.jl * https://github.com/under-Peter/OMEinsum.jl * https://github.com/mcabbott/TensorCast.jl

That seems strange, if I'm reading it correctly.

In the author's packages, integers which label dimensions etc. have gradient `nothing`, but arrays which happen to contain integers do not signal anything:

    julia> Zygote.gradient((x,d)->abs(sum(sin,fft(x,d))), [2 2; 0 0], 2)
    (Complex{Float64}[-0.346356+0.0im 1.65364+0.0im; -2.0+0.0im 0.0+0.0im], nothing)

A tensor product of two matrices is just this, with 4 loose ends:

    ---A---
    ---B---

A_ab B_cd = (A⊗B)_acbd. If a,c=1...N then you can think of (ac) as being an index with N^2 possible values, which is what the double line now means... if you like, you may compact it a bit:

    ===(A⊗B)===

The thing about RPN is that it worked well as a sort-of "spoken not written" language. It's a clear and usable way to give instructions one-at-a-time to a calculator with a 1-line display.

But written out on paper, it's very hard to see what's going on at a glance, which is why nobody does this on paper.

Almost anyone programming in mathematical expressions of reasonable complexity today will spend much more time reading (and triple-checking!) these expressions than typing them. That's a strong reason to like infix notation, and in general notation close to what you would use on paper.

Maybe I phrased that poorly: you can indeed write the loops yourself (without a speed penalty) as advertised, to do things which are hard/ugly/impossible in numpy. But very often you can do this without dropping all the way down to explicit pointer arithmetic, where you'd have to care about base-0/1 and off-by-one errors (again without penalty). This in-between zone is what I was trying to suggest is worth exploring. See for example https://julialang.org/blog/2016/02/iteration for a taste.

Someone more polyglot than I am will have to comment on C++ comparisons.