Low rank factorized AbstractMatrix

General Usage Performance
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).

1 Like
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 Likes
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)
1 Like
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 Interfaces · The Julia Language. 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.

2 Likes
7

Okay. This provides more discussion: Interfaces for Abstract Types · Issue #6975 · JuliaLang/julia · GitHub

1 Like
8

You can use LowRankMatrix in LowRankApprox.jl

4 Likes
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 _.

1 Like
11

I’ll give it a try.

12

GitHub - timholy/WoodburyMatrices.jl: Support for the Woodbury matrix identity for Julia should also support that.

1 Like
13

Sorry to revive an old topic. But I found myself needing this. I was wondering if there are updated recommendations to achieve this?

Also, would it be a good idea to break LowRankMatrix apart into its own package?

Thanks!

14

Would it be okay if I did this, with proper attribution to the origin of the code? With some minor changes, probably.

15

Yes, it’ll be ok, although what’s the point, instead of contributing changes to the package? It’s MIT expat licensed, so one may freely copy with attributions.

1 Like
16

LowRankApprox would be overkill for me. I just need a simple LowRankMatrix type.

To avoid code duplication, the possibility I would suggest is that LowRankApprox depends on a putative LowRankMatrix package.

17

Consider also if you want to change the name. It may be better called MatrixProduct instead of LowRankMatrix since the data structure is not specific for low-rank matrices, although the product of a tall skinny matrix and a fat short matrix is a common use case.

18

True, but the idea is that the implementation of things like matrix product exploit the tall / skinny shapes of a low-rank product. So I think LowRankMatrix makes sense.

2 Likes
19

Just use LazyArrays.jl which has ApplyArray(*, A, B)