In Julia there are some subtypes of AbstractMatrix, which are used as "lazy wrapper" objects for objects of type AbstractMatrix. These inlcude Transpose, Adjoint, Symmetric, Hermitian, UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular and maybe others.
The purpose of those objects is primarily to avoid unnecessary moving or reorganizing element data of the underlying matrices and instead give a different view to the parent data. The use of the wrappers can improve performance, if there are specialized algorithms, which implicitly make use of the reorganized data, as is the case for the BLAS library or some methods for sparse matrices.
It is possible to construct arbitrary nestings of such wrapped objects. Two nested wrapped matrices may be different in a strict sense - i.e. they have different types - while the mathematical meaning is identical.
For example we have Transpose(UpperTriangular(A)) == LowerTriangular(Transpose(A)).
Unfortunatly deeper nested wrapped objects are not well supported by specialized algorithms, presumably because there is an undefinite number of different nested types. At the same time the existence of simplifications shows, that the required number of specialized algorithms is smaller. If we have a mathematical function f(::Transpose{T,UpperTriangular{T,S}}), we don't need an extra f(::LowerTriangular{T,Transpose{T,S}}) because the mathematical meaning of the arguments is identical.
The purpose of this writing is to show, that there is only a finite number of mathematical meanings for objects, which have been constructed by application of a arbitrary number of "lazy wrappers" to any square matrix. The approach is constructive by presenting a "UniversalWrapper", a parameterized subtype of AbstractMatrix, which covers all the cases of combinations of wrappers. The number of parameter combinations is 18 if restricted to real matrices and 75 in general. Some of the potential parameter combinations are unused, so the valid group size is less.
We must admit, that not all kinds of "lazy wrappers" are supported by this approach. For example view respectively SubArray can't be handled; maybe with the exception of view(A, I, I) with identical row and column indices. I think, an unrelated approach for SubArray is required. Also "structured types" Diagonal, BiDiagonal, Tridiagonal, ... are not covered.
- Equivalence Relation:
A == B <=> size(A) == size(B) & ∀ i,j ( A[i,j] === B[i,j] || A[i,j] == B[i,j]) - Well-defined:
f respects equivalence relation <=> A == B => f(A) == f(B)
The === in the definition of the equivalence relation is due to the existence of NaN, which prevents == from being an equivalence relation on Float64/32. For matrices not containing NaN this definition coincides with the definition of Juliamethod ==(A::AbstractMatrix, B::AbstractMatrix). In the rest of the text we ignore the existence of NaN.
Each "lazy wrapper" can be viewed as a transformation of a subset of the space of matrices into itself. We ignore the algebraic structure of this space of matrices here (that means addition and multiplication). Actually a matrix is reduced to a mapping of integer indices [i,j] to numbers. The numbers may be real or complex.
Such sets of transformations with composition as binary operation form a Monoid, a special kind of Endomorphism, which is close to a group. While a unity element exists (the identity transformation), inverses don't exist for all transformations.
The size of the monoid depends on the subset of matrices, on which the transformations operate. If we restrict to real matrices, some of the transformations for complex matrices are not relevant. The same is true, if we restrict to complex matrices with real diagonal, because then the transformations, which modify diagonal elements (like Hermitian) are redundant. We can also reduce the size, if we ignore the transformations transpose and symmetric, which may be mathematically uninteresting in the complex case.
In Julia a difference between "functions" and "constructors" exists. Constructors share the name with an abstract or concrete type and there is an unwritten law, that requires Constructor(args...) isa Constructor. Unfortunately this implies, that nested application of wrapper constructors cannot apply obvious simplifications. For example Transpose(Transpose(A)) !== A. Therfore wrapper functions have been introduced, which fall back to corresponding constuctors in simple cases, but which perform obvious simplifications in the nested case. For example for A isa Matrix while transpose(A) === Transpose(A) we have transpose(transpose(A)) === A.
There is the idiosyncrasy of Symmetric/Hermitian, where the constructurs and corresponding wrapper functions require the extra argument uplo in [:U, :L]. We introduce separate wrapper functions for the both cases.
Because we want to enforce simplifications in a consistent way, we introduce "wrapper functions" for all of the supported types. The following list shows the relationships:
| Wrapper Type | Constructor(A) | Wrapper Function |
|---|---|---|
| Transpose | Transpose(A) | transpose |
| Adjoint | Adjoint(A) | adjoint |
| Symmetric | Symmetric(A, :U) | uppersymm |
| Symmetric(A, :L) | lowersymm | |
| Hermitian | Hermitian(A, :U) | upperherm |
| Hermitian(A, :L) | lowerherm | |
| UpperTriangular | UpperTriangular(A) | uppertria |
| LowerTriangular | LowerTriangular(A) | lowertria |
| UnitUpperTriangular | UnitUpperTriangular(A) | unituppertria |
| UnitLowerTriangular | UnitLowerTriangular(A) | unitlowertria |
| diagonal | ||
| conjugate | ||
| identity | ||
| one |
Most of the wrapper functions are new. For the last two lines there exists no wrapper types yet. All wrapper types shall be covered by a new universal type, with conversions from and partially to the existing types.
==is an equivalence relation on the set of square matrices not containingNaNelements- all wrapper functions are well-defined on this set
Let n a positive integer and A an abstract matrix with A[i,j] isa Number and !isnan(A[i,j]) for all i, j in 1:n.
Let up, lo in -2:2 and di in 0:2.
Define U{up,di,lo}(A) as an abstract matrix with
U{up,di,lo}(A)[i,j] :=
if i < j ?
up == 1 -> A[i,j]
up == 2 -> A[j,i]
up == 0 -> 0
up == -1 -> conj(A[i,j])
up == -2 -> conj(A[j,i])
elseif i > j ?
lo == 1 -> A[i,j]
lo == 2 -> A[j,i]
lo == 0 -> 0
lo == -1 -> conj(A[i,j])
lo == -2 -> conj(A[j,i])
else # i == j
di == 0 -> 1
di == 1 -> A[i,i]
di == -1 -> conj(A[i,i])
di == 2 -> real(A[i,i])
end
If we restrict to real element types, which have conj(a) == real(a) == a, all negative parameters and di == 2 are not required.
In this case the defining function is simplified correspondingly.
With the definition of U the standard wrapper functions correspond to certain transformations in the following way:
| Wrapper Type | Constructor(A) | Wrapper Funct | Transformation |
|---|---|---|---|
| Transpose | Transpose(A) | transpose | U{2,1,2} |
| Adjoint | Adjoint(A) | adjoint | U{-2,-1,2} |
| Symmetric | Symmetric(A, :U) | uppersymm | U{1,1,2} |
| Symmetric(A, :L) | lowersymm | U{2,1,1} | |
| Hermitian | Hermitian(A, :U) | upperherm | U{1,2,-2} |
| Hermitian(A, :L) | lowerherm | U{-2,2,1} | |
| UpperTriangular | UpperTriangular(A) | uppertria | U{1,1,0} |
| LowerTriangular | LowerTriangular(A) | lowertria | U{0,1,1} |
| UnitUpperTriangular | InitUpperTriangular(A) | unituppertria | U{1,0,0} |
| UnitLowerTriangular | UnitLowerTriangular(A) | unitlowertria | U{0,0,0} |
| diagonal | U{0,1,0} | ||
| conjugate | U{-1,-1,-1} | ||
| identity | U{1,1,1} | ||
| one(A) | U{0,0,0} |
Composing two transformations U1 = U{up1,di1,lo1} and U2 = U{up2,di2,lo2} is defined by
(U1 ∘ U2)(A) := U1(U2(A))
There is always a U3 = U{up3,di3,lo3} with U3 == U1 ∘ U2.
The type-defining integers up3,di3,lo3 may be calculated using the scheme defined in implementation.
Let na be a positive integer, 𝕂 = ℝ or 𝕂 = ℂ.
The previously defined set of real matrix transformations of form a finite submonoid of the monoid of all endomorphisms of 𝕂^n,n. The unit element is U{1,1,1}. The unique zero element is U{0,0,0}.
The monoid is generated by {unituppertria, uppertria, unitlowertria, lowertria, transpose, uppersymmetric, lowersymmetric} in the case of 𝕂 = ℝ and {unituppertria, uppertria, unitlowertria, lowertria, transpose, adjoint, uppersymmetric, lowersymmetric, upperherm, lowerherm} in the case of 𝕂 = ℂ.