Sketch of linear relations in julia

linrel.jl

#A*x = 0
using LinearAlgebra
struct HRep{T}
   n :: Int64 #input size
   A :: AbstractArray{T}
end
struct VRep{T}
   n :: Int64 #input size
   A :: AbstractArray{T}
end

function vrep(x::HRep) 
    VRep(x.n, nullspace(x.A))
end
function hrep(x::VRep)
    HRep(x.n, nullspace(x.A')')
end


function meet(x::HRep, y::HRep)
    @assert x.n == y.n
    HRep(x.n, [x.A ; y.A])
end
function meet(x::VRep, y::VRep)
    @assert x.n == y.n
    xh = hrep(x)
    yh = hrep(y)
    HRep(x.n, [xh.A ; yh.A])
end
function meet(x::HRep, y::VRep)
    @assert x.n == y.n
    yh = hrep(y)
    HRep(x.n, [x.A ; yh.A])
end
function meet(x::VRep, y::HRep)
    @assert x.n == y.n
    xh = hrep(x)
    HRep(x.n, [xh.A ; y.A])
end


function join(x::VRep, y::VRep)
    @assert x.n == y.n
    VRep(x.n, [x.A y.A])
end
function join(x::HRep, y::HRep)
    yv = vrep(y)
    join(x,yv)
end
function join(x::VRep, y::HRep)
        xv = vrep(x)
        join(xv,y)
end
function join(x::HRep, y::HRep)
        xv = vrep(x)
        join(xv,y)
end
    
function rid(n)
    VRep(n, [ Matrix(I,n,n) ; Matrix(I,n,n) ]  )
end


function rsub(r :: VRep,p :: HRep)
    all(isapprox.(p.A * r.A , 0; atol=eps(Float64), rtol=0))
end

function rsub(r :: VRep, p :: VRep)
    p = hrep(p)
    all(isapprox.(p.A  * r.A , 0; atol=eps(Float64), rtol=0))
end

function heq(p,q)
   rsub(p,q) && rsub(q,p) 
end


function converse(p :: VRep)
   n, g = size(p.A)
    VRep( n - p.n  , [p.A[p.n + 1 : end, :] ;  p.A[1 : p.n, :]  ]) 
end

function converse(p :: HRep)
   c, n = size(p.A)
   HRep( n - p.n  , [p.A[:, p.n + 1 : end]  p.A[:, 1 : p.n]]) 
end

function top(n)
   VRep(n, Matrix(I, n, n)) 
end 
function bottom(n)
    HRep(n, Matrix(I, n, n)) 
end


function compose(x::HRep, y::HRep)
    cx, nx = size(x.A)
    cy, ny = size(y.A)    
    na = x.n
    nb = nx - x.n # which should equal y.n
    nc = ny - y.n
    B = [                x.A          zeros( cx, ny - y.n) ; 
         zeros( cy, x.n  )      y.A                        ]
    C = nullspace(B)
    return VRep([  C[1 : n.x, :]          ; 
                   C[ nx+y.n + 1:end, :] ])
end

I'm a little confused as to how compose(::HRep, ::Hrep) computes the composition of two relations. Could you explain that a little more? In particular, I don't see where the codomain part of x interacts with the domain part of y.

Well, I may not have implemented it correctly.

The h representation is a vertical stack of vectors for which the subspace you're describing dots to 0. There are 3 spaces in discussion, one of size na, one of size nb, and one of size nc.

x is cx tall stack of vectors na + nb) long, y is a stack of cy vectors (nb + nc) long. I extend x and y with zeros so that they are both (na + nb + nc) long and then vertically stack those together to create the matrix B which is size (cx + cy) x (na + nb + nc) .

I convert it to Vrep by taking the nullspace. These are the generators of the subspace of interest. I then project out the interior nb variables, which is the definition of relational composition (merely that values exist).

I tried doing something a little different from my blog post http://www.philipzucker.com/linear-relation-algebra-of-circuits-with-hmatrix/ in that I'm using linear relations here rather than affine. I think it will come out cleaner, and then one can embed the affine version in this by a standard trick of extending all vectors by 1 dimension that corresponds to the constant offset part. See homogenous coordinates for more https://en.wikipedia.org/wiki/Homogeneous_coordinates

Oh yeah. I DEFINITELY didn't implement it correctly

I tried to fix it up a little. I still stand by the text of my reasoning above, but it may not be correct by not actually being julia syntax and being off on indexing calculation.

I tried formatting that matrix a little better so you can see where the overlap is. That B matrix is not cut like standard block matrix notation. The x.A and y.A blocks overlap for nb horizontal entries

I think I worked it out on paper just now,

f: R^n -> R^m
g: R^m -> R^k

x,y in f iff Ax = By
y,z in g iff Cy = Dz

[A B] [x]  = 0   & [C D] [y]  = 0 
      [y]                [z] 


[ A B 0 ]  [x]  = [0]
[ 0 C D ]  [y]    [0]
           [z]

So you want the nullspace([A B 0; 0 C D]) and then project out the dimensions of dim(A)+1:dim(A)+1+dim(B) where dim(A) is the number of columns.

Yup. We're on the same page. I think what you wrote is much clearer giving the blocks explicit names

Is your B matrix the same as my

[ A B 0 ]
[ 0 C D ]

?

Yes
Edit: modulo me screwing it up

Asynchronous Chat! I'm not seeing your posts until I submit mine!

Truly this is how github was meant to be used

Cool, thanks for helping me understand it. I wish LinRel was taught in Linear Algebra Class, I had no idea about it until a few months ago.

Me either basically. Systems of linear equations / subspaces is not novel certainly. That there is a high level algebra for it is pretty interesting. My impression of the category theorist version of this is that they want to make a synthetic version of this (without numbers), which is very strange to my sensibilities.

I think one of the cool accidental upsides of the ACT approach to linear algebra is that it naturally captures the informal "block matrix notation" that everyone uses but just kinda wills into existence without proving all those boring foundations details.

Is block notation informal?

It is often taught informally. Unlike scalar notation where the professor spends multiple lectures building the equivalence between matrices and linear operators. Block notation is introduced as “it’s like regular notation but the entries are matrices instead of numbers”

Is there a range version of this like:

f: R^n -> R^m
g: R^m -> R^k

x,y in f iff x = Au and y=Bu
y,z in g iff y = Cv and z = Dv

[ A  0 ]  [u]  = [x]
[ B -C ]  [v]    [0]
[ 0  D ]         [z]

And then projecting out the middle rows?

This does not seem that familiar to me.
My method for composing range versions (which I called VReps) was to first convert them to the hrep

Are you requiring y to be 0? I'm not sure why one would. And if you are, then there is an implicit solve there
Does what you're suggesting have the identity relation as an identity?

Projecting out the middle rows would be just stacking A and D?

function compose(x::VRep, y::VRep)
   x1 = hrep(x)
   y1 = hrep(y)
   compose(x1,y1)
end