Low rank factorized AbstractMatrix


#1

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:

B=sparse(u)*sparse(u')

The sparse approach is better but not optimal

julia> x=randn(n1+n2);
julia> @btime A*x;
  345.225 μs (1 allocation: 7.94 KiB)
julia> @btime B*x;
  9.531 μs (1 allocation: 7.94 KiB)
julia> @btime u*(u'*x);
  1.044 μs (3 allocations: 7.97 KiB)

and it uses that the vector is sparse (which is not always the case in the application).


#2

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

#3

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:

#4

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 you can define it as:

Base.getindex(m::LowRankMatrix, i::Int, j::Int) = sum(x -> x[i] * x[j], m.vecs)

#5

Thanks a lot.

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.


#6

Interfaces are by convention in Julia, not a language feature (yet). You can take a look at examples and copy them. Here is another example: https://github.com/JuliaArrays/OffsetArrays.jl/blob/master/src/OffsetArrays.jl.

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.


#7

Okay. This provides more discussion: https://github.com/JuliaLang/julia/issues/6975


#8

You can use LowRankMatrix in LowRankApprox.jl


#9

Indeed, this works very nicely:

julia> using LowRankApprox
julia> C=LowRankApprox._LowRankMatrix(u,u);
julia> @btime C*x;
  1.526 μs (2 allocations: 8.03 KiB)

Since you are also the package author: Is there no exported creator such I can specify the factorization directly?


#10

I don’t think so, but I think the use of _ here was a mistake. I’ll happily merge a PR that removes the _.


#11

I’ll give it a try.


#12

https://github.com/timholy/WoodburyMatrices.jl should also support that.