Operator example: trace

This page illustrates the "linear operator" feature of the Julia package LinearMapsAA.

This page comes from a single Julia file: 02-trace.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: 02-trace.ipynb, or open it in binder here: 02-trace.ipynb.

Setup

Packages needed here.

using LinearMapsAA
using LinearAlgebra: tr, I
using InteractiveUtils: versioninfo

Overview

The "operator" aspect of this package may seem unfamiliar to some readers who are used to thinking in terms of matrices and vectors, so this page describes a simple example: the matrix trace operation. The trace of a $N × N$ matrix is the sum of its $N$ diagonal elements. We tend to think of this a function, and indeed it is the tr function in the LinearAlgebra package. But it is a linear function so we can represent it as a linear operator $𝒜$ that maps a $N × N$ matrix into its trace. In other words, $𝒜 : \mathbb{C}^{N × N} \mapsto \mathbb{C}$ is defined by $𝒜 X = \mathrm{tr}(X)$. Note that the product $𝒜 X$ is not a "matrix vector" product; it is a linear operator acting on the matrix $X$.

(Note that we use a fancy unicode character $𝒜$ here just as a reminder that it is an operator; in practical code we usually just use A.)

The LinearMapsAA package can represent such an operator easily. Here is the definition for $N = 5$:

N = 5
forw(X) = [tr(X)] # forward mapping function
𝒜 = LinearMapAA(forw, (1, N*N); idim = (N,N), odim = (1,))

LinearMapAO: 1 × 25 odim=(1,) idim=(5, 5) T=Float32 Do=1 Di=2
map = 1×25 LinearMaps.FunctionMap{Float32,false}(#16; issymmetric=false, ishermitian=false, isposdef=false)

The idim argument specifies that the input is a matrix of size N × N and the odim argument specifies that the output is vector of size (1,).

One subtlety with this particular didactic example is that the ordinary trace yields a scalar, but LinearMaps.jl is (understandably) designed exclusively for mapping vectors to vectors, so we use [tr(X)] above so that the output is a 1-element Vector. This behavior is consistent with what happens when one multiplies a 1 × N matrix with a vector in $\mathbb{C}^N$.

Here is a verification that applying this operator to a matrix produces the correct result:

X = ones(5)*(1:5)'
𝒜 * X, tr(X), (N*(N+1))÷2

([15.0], 15.0, 15)

Although $𝒜$ here is not a matrix, we can convert it to a matrix (at least when $N$ is sufficiently small) to perhaps understand it better:

A = Matrix(𝒜)
A = Int8.(A) # just for nicer display

1×25 Matrix{Int8}:
 1  0  0  0  0  0  1  0  0  0  0  0  1  0  0  0  0  0  1  0  0  0  0  0  1

The pattern of 0 and 1 elements is more obvious when reshaped:

reshape(A, N, N)

5×5 Matrix{Int8}:
 1  0  0  0  0
 0  1  0  0  0
 0  0  1  0  0
 0  0  0  1  0
 0  0  0  0  1

Adjoint

Although this is largely a didactic example, there are optimization problems with trace constraints of the form $𝒜 X = b$. To solve such problems, often one would also need the adjoint of the operator $𝒜$.

Mathematically, and adjoint is a generalization of the (Hermitian) transpose of a matrix. For a (bounded) linear mapping $𝒜$ between inner product space $𝒳$ with inner product $⟨ \cdot, \cdot \rangle_𝒳$ and inner product space $𝒴$ with inner product $⟨ \cdot, \cdot \rangle_𝒴,$ the adjoint of $𝒜$, denoted $𝒜'$, is the unique bound linear mapping that satisfies $⟨ 𝒜 x, y \rangle_𝒴 = ⟨ x, 𝒜' y \rangle_𝒳$ for all $x ∈ 𝒳$ and $y ∈ 𝒴$.

Here, let $𝒳$ denote the vector space of $N × N$ matrices with the Frobenius inner product for matrices: $⟨ A, B \rangle_𝒳 = \mathrm{tr}(A'B)$. Let $𝒴$ simply be $\mathbb{C}^1$ with the usual inner product $⟨ x, y \rangle_𝒴 = x_1^* y_1$.

With those definitions, one can verify that the adjoint of $𝒜$ is the mapping $𝒜' c = c_1 \mathbf{I}_N$, for $c ∈ \mathbb{C}^1$, where $\mathbf{I}_N$ denotes the $N × N$ identity matrix.

Here is the LinearMapAO that includes the adjoint:

back(y) = y[1] * I(N) # remember `y` is a 1-vector
𝒜 = LinearMapAA(forw, back, (1, N*N); idim = (N,N), odim = (1,))

LinearMapAO: 1 × 25 odim=(1,) idim=(5, 5) T=Float32 Do=1 Di=2
map = 1×25 LinearMaps.FunctionMap{Float32,false}(#6, #8; issymmetric=false, ishermitian=false, isposdef=false)

Here is a verification that the adjoint is correct (very important!):

@assert Matrix(𝒜)' == Matrix(𝒜')
Int8.(Matrix(𝒜'))

25×1 Matrix{Int8}:
 1
 0
 0
 0
 0
 0
 1
 0
 0
 0
 ⋮
 0
 0
 1
 0
 0
 0
 0
 0
 1

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

12-element Vector{SubString{String}}:
 "Julia Version 1.10.3"
 "Commit 0b4590a5507 (2024-04-30 10:59 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()

Status `~/work/LinearMapsAA.jl/LinearMapsAA.jl/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [e30172f5] Documenter v1.4.1
  [7a1cc6ca] FFTW v1.8.0
  [71a99df6] ImagePhantoms v0.8.0
  [7031d0ef] LazyGrids v1.0.0
  [599c1a8e] LinearMapsAA v0.12.0 `~/work/LinearMapsAA.jl/LinearMapsAA.jl`
  [98b081ad] Literate v2.18.0
  [170b2178] MIRTjim v0.24.0
  [b77e0a4c] InteractiveUtils
  [37e2e46d] LinearAlgebra

This page was generated using Literate.jl.