Here I show an extreme example of matrix-free calculations done completely in the abstract using Ocaml modules and functors. Previously I showed how we can use Ocaml’s module capabilities to abstract away matrices from numerical linear algebra. This enabled calculation of interesting spectral values like operator adjoints and operator norms without once refering to the operator’s “matrix.” I gave an explicit example of polynomial differentiation which could be represented as a matrix, but instead I represented it as recursion over a syntax tree. Here I will define a vector space which is infinite dimensional and define an operator on this vector space which can not be represented as a matrix, and then I compute its operator norm using exactly the same code used in the polynomial case.
I won’t go into great detail of how the module components work; I have done that elsewhere. I will just state briefly how I defined the critical operations on the vector space in ocaml, and follow with results and verification.
I want this vector space to encapsulate “the space of smooth functions.” Thus I defined the vector type as follows:
(*The "vect" type is the type of functions
with one floating point argument and floating point output*)
type vect = (float -> float)
The basic axioms of vector spaces only stipulate that one vector need exist, the null vector. For functions that will be the constant function zero:
(* nullvector : vect *)
let nullvector = fun x -> 0.0
I also added the requirement that at least one basis be provided to the modules, otherwise virtually nothing can be done in a truly computable fashion.
let rec poly i x = if i==0 then (1.0) else ( x *. (poly (i-1) x))
(* basis i : vect *)
let basis i = poly i
The basis here can be virtually anything though, the only true requirement that they are all of type vect. I chose polynomials here for the approximation theorems and the simplicity.
A vector space is a set of vectors, a field, and operations defined between the vectors and field elements. Here I show how I defined vector addition, scalar multiplication, and the inner product.
(* Addition of functions: (f1+f2)(x) = f1(x) + f2(x) *)
let add f1 f2 = fun x -> (f1 x) +. (f2 x)
(* Scalar multiplication: (s*f)(x) = s*f(x) *)
let scalarmul s f = fun x -> s *. (f x)
(* Inner product of two functions f1,f2:
integral of f1*f2 on interval [-1,1] *)
let innerprod f1 f2 =
let f = fun x -> (f1 x) *. (f2 x) in
Quad.gkint f a b tol
let norm p = (sqrt (innerprod p p))
Note that the Quad.gkint above is an adaptive quadrature routine I wrote myself. You may find it here.
With the functions in hand we just need to put it into a module and input that module into the functors which defines the matrix-free computations. I have done that here and I have called the module containing the above functions FunctionHilbert If the module defining field operations is called “Real” then the functor may be called as below:
(* Define the module for computations on operators of type vect -> vect *)
(* Note the "Real" module input here is just a module that defines the field
operations for floating point numbers. *)
module MyOp = MakeOperatorSpace (Real) (FunctionHilbert) (FunctionHilbert)
Now that I have a vector space I can define a linear operator of type vect->vect. I have picked a simple operator on the space of functions:
(* "a" modulates "f" by the trig function "sin" *)
let a f = fun x -> (sin x) *. (f x)
Note that this linear operator takes an input function “f” and outputs another function which is “f” modulated by the trig function “sin.” This has no easy matrix representation the way polynomial derivative operator does. I chose this example just to illustrate that the operator norm code works even in this extreme case.
First we can use the MyOp module to compute the operator norm of a. I have it run 150 iterations and on 10 basis functions, which was the point it seemed to converge to a consistent number. (Note: We know this will converge because of approximation theorems concerning polynomials: they are dense in the set of smooth functions, hence my deliberate selection of polynomials as basis functions)
let opnorm = MyOp.opnorm 150 10 a
let () = print_endline (Real.to_string opnorm)
with corresponding output.
>0.830681896978
The amazing thing is that this uses exactly the same code as my polynomial derivative example. I changed nothing to achieve this much more general and interesting computation
The above operator norm calculation was done by appealing to linear algebraic properties, specifically by computing the spectral radius of an operator. We can verify that the above number is accurate by appealing directly to the definition of operator norm, which is posed as an optimization problem. Rather than compute this by hand, I just wrote a simple Julia script to solve the optimization problem:
using Polynomials;
using Optim;
function opnorm(a)
p=Poly(a);
s(x) = sin(x)*sin(x);
q(x) = p(x)*p(x);
f(x) = s(x)*q(x);
#Use Julia's built-in adaptive quadrature.
num=quadgk(f,-1.0,1.0)[1];
den=quadgk(q,-1.0,1.0)[1];
return -sqrt(num/den);
end
println( abs(optimize(opnorm,ones(10)).f_minimum ));
which outputs
>0.8272431868114323
which is virtually the same value computed for the operator norm, just directly from its definition rather than by appealing to linear algebra.
The linear operator I defined here is known as a “multiplication operator.” Furthermore this operator is symmetric (equals its adjoint). Thus the operator norm here coincides with the spectral radius of the operator itself. Thus appealing to Paul Halmos in Introduction to Hilbert Space and the Theory of Spectral Multiplicity I know that the entire spectrum of the operator will just be the entire range of sin. The spectral radius of this is just the maximum of the spectrum, which in this case is
sin(1) = 0.8414709848...
What amazed me in working through this example was how simple it was to write such abstract code in Ocaml. This vector space is an especially badly behaved example in math because it is infinite dimensional and many linear operators defined on it are unbounded, and so have no operator norm (or, depending on your viewpoint, have infinite operator norm). Furthermore I constructed specifically an operator with no matrix representation and successfully used exactly the same code to compute its operator norm as what I used to compute the operator norm of the polynomial derivative operator here.
There are a few caveats to point out. The inner product I defined is not a true inner product because it works by pointwise sampling functions in order to compute their integrals. Therefore we could “trick” the inner product into thinking a function is constantly 0 even when it is not, just by explicitly defining the function to be 0 at the sampling points of the integration routine. This could be overcome theoretically by instead considering the type of vectors to be an equivalence class of functions, but this has questionable value since it makes the code more complex to overcome a very unlikely scenario.
Another issue is the way I have use lambdas. Since the vector arithmetic is returning a lambda, we could cause stack overflows by doing too much arithmetic - say in a deeply nested loop. This could be overcome by defining a custom function type with permissible operations, that way allowing periodic simplification to keep in check the complexity of deeply nested arithmetic.