Basic question: Is there an AbstractMatrix which can represent a low-rank factorized matrix?

Example: I have the rank-one matrix:

n1=900; n2=100;
u=[zeros(n1);randn(n2)];
A=u*u';

I want to use A in functions which take AbstractMatrix as input, and use many matrix vector products (but also sum matrices together). Is there a way to represent u*u', without multiplying it out?

Current solution: I use that u has a sparsity pattern and represent A as a sparse matrix:

Oh well, the beauty of Julia is that you can code such a matrix up in a few minutes and get great performance! Here is an example. You can improve on the mul! and add more functionality if you want but this should get you started:

struct LowRankMatrix{T, Tvecs <: Tuple{Vararg{AbstractVector{T}}}} <: AbstractMatrix{T}
vecs::Tvecs
end
function Base.size(m::LowRankMatrix, i)
(i == 1 || i == 2) && return length(m.vecs[1])
return 1
end
function Base.size(m::LowRankMatrix)
return (size(m, 1), size(m, 2))
end
LowRankMatrix(v::AbstractVector) = LowRankMatrix((v,))
Base.:+(m1::LowRankMatrix, m2::LowRankMatrix) = LowRankMatrix((m1.vecs..., m2.vecs...))
using LinearAlgebra
function LinearAlgebra.mul!(c::AbstractVector, m::LowRankMatrix, b::AbstractVector)
c .= 0
for vec in m.vecs
axpy!(dot(vec, b), vec, c)
end
return c
end

Usage

julia> v = [1.0, 2.0];
julia> m1 = LowRankMatrix(v);
julia> m2 = LowRankMatrix(v);
julia> b = rand(2);
julia> m1 * b == dot(v, b) * v
true
julia> (m1 + m2) * b == m1 * b + m2 * b == 2 * (m1 * b)
true

Thanks. How do I know which methods I need to reimplement to satisfy what is expected from an AbstractMatrix? It seems like it could be a lot. Can I make it use some fallback for methods I have not implemented?

For example: Your nice type throws error when accessing individual elements:

julia> m1[1,2]
ERROR: getindex not defined for LowRankMatrix{Float64,Tuple{Array{Float64,1}}}
Stacktrace:

How do I know which methods I need to reimplement to satisfy what is expected from an AbstractMatrix ?

You may not need all the functions of AbstractMatrix. Implementing them may be a cool project, but could also be overkill. You can take a look at https://github.com/JuliaArrays/FillArrays.jl/tree/master/src for the list of functions you can choose to implement.

m1[1,2] is actually getindex with 2 arguments.

julia> @which m1[1,2]
getindex(A::AbstractArray, I...) in Base at abstractarray.jl:903

So there are no "standardized methods associated with AbstractMatrix"? (Like an interface in Java, I guess.) Your solution would solve the problem for now but does not seems so reliable from a system design point of view. I want to hand over my LowRankMatrix to some code I do not have control over. If I now reverse engineer and try to figure out which AbstractMatrix-functions that code uses, it would work now for the cases I test, but it does not rule out that they change the code to use more functions in the future, or that I use other parameters in their code which lead to other functions.

The reference you give is to an example. I was hoping for a standard. Anyway thanks for the help.

Edit: The docs actually document a number of functions typically defined for arrays https://docs.julialang.org/en/v1/manual/interfaces/index.html#man-interface-array-1. And not every array type actually defines all of these. For example immutable arrays in StaticArrays.jl don’t have setindex!, or actually they do but it errors. There is also the question of how much linear algebra you want to support.