Installation

Add ParametricOperators.jl as a dependency to your environment.

To add, either do:

julia> ]
(v1.9) add ParametricOperators

OR

julia> using Pkg
julia> Pkg.activate("path/to/your/environment")
julia> Pkg.add("ParametricOperators")

Simple Operator

Make sure to include the package in your environment

using ParametricOperators

Lets start by defining a Matrix Operator of size 10x10:

A = ParMatrix(10, 10)

Now, we parametrize our operator with some weights θ:

θ = init(A)

We can now apply our operator on some random input:

x = rand(10)
A(θ) * x

Gradient Computation

Make sure to include an AD package in your environment

using Zygote

Using the example above, one can find the gradient of the weights or your input w.r.t to some objective using a standard AD package:

# Gradient w.r.t weights
θ′ = gradient(θ -> sum(A(θ) * x), θ)

# Gradient w.r.t input
x′ = gradient(x -> sum(A(θ) * x), x)

Chaining Operators

We can chain several operators together through multiple ways

Compose Operator

Consider two matrices:

L = ParMatrix(10, 4)
R = ParMatrix(10, 4)

We can now chain and parametrize them by:

C = L * R'
θ = init(C)

This allows us to perform several operations such as

using Zygote
using LinearAlgebra

x = rand(10)
C(θ) * x
gradient(θ -> norm(C(θ) * x), θ)

without ever constructing the full matrix LR', a method more popularly referred to as LR-decomposition.

Kronecker Operator

Kronecker Product is a most commonly used to represent the outer product on 2 matrices.

We can use this to describe linearly separable transforms that act along different dimensions on a given input tensor.

For example, consider the following tensor:

T = Float32
x = rand(T, 10, 20, 30)

We now define the transformation that would act along each dimension. In this case, a Fourier Transform.

Fourier Transform Example

Fx = ParDFT(T, 10)
Fy = ParDFT(Complex{T}, 20)
Fz = ParDFT(Complex{T}, 30)

We can now chain them together using a Kronecker Product:

F = Fz ⊗ Fy ⊗ Fx

Now, we can compute this action on our input by simply doing:

F * vec(x)

    Sx = ParIdentity(T, 10)
    Sy = ParMatrix(T, 20, 20)
    Sz = ParMatrix(T, 30, 30)

    S = Sz ⊗ Sy ⊗Sx
    θ = init(S)

    S(θ) * vec(x)