# Low rank factorized AbstractMatrix

**URL:** <https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838>\
**Category:** Performance\
**Created:** [January 20, 2019, 11:03am UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838 "2019-01-20T11:03:30Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [January 20, 2019, 11:03am UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/1 "2019-01-20T11:03:30Z")

</div>

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

Example: I have the rank-one matrix:

```julia
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:

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

```

The sparse approach is better but not optimal

```julia
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).

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [January 20, 2019, 12:27pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/2 "2019-01-20T12:27:08Z")

</div>

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:

```julia
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
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

```

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [January 20, 2019, 12:45pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/3 "2019-01-20T12:45:23Z")

</div>

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
julia> m1[1,2]
ERROR: getindex not defined for LowRankMatrix{Float64,Tuple{Array{Float64,1}}}
Stacktrace:

```

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [January 20, 2019, 12:57pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/4 "2019-01-20T12:57:07Z")

</div>

> 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](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
julia> @which m1[1,2]
getindex(A::AbstractArray, I...) in Base at abstractarray.jl:903

```

So you can define it as:

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

```

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [January 20, 2019, 1:06pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/5 "2019-01-20T13:06:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [January 20, 2019, 1:10pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/6 "2019-01-20T13:10:08Z")

</div>

> [@jarl](#):
>
> So there are no "standardized methods associated with `AbstractMatrix` "? (Like an interface in Java, I guess.)

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](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](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.

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [January 20, 2019, 1:19pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/7 "2019-01-20T13:19:23Z")

</div>

Okay. This provides more discussion: [Interfaces for Abstract Types · Issue #6975 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/6975)

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 20, 2019, 3:33pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/8 "2019-01-20T15:33:15Z")

</div>

You can use `LowRankMatrix` in LowRankApprox.jl

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [January 20, 2019, 5:00pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/9 "2019-01-20T17:00:56Z")

</div>

Indeed, this works very nicely:

```julia
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?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 20, 2019, 5:03pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/10 "2019-01-20T17:03:32Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [January 20, 2019, 5:15pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/11 "2019-01-20T17:15:35Z")

</div>

I’ll give it a try.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [January 20, 2019, 5:37pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/12 "2019-01-20T17:37:40Z")

</div>

[GitHub - timholy/WoodburyMatrices.jl: Support for the Woodbury matrix identity for Julia](https://github.com/timholy/WoodburyMatrices.jl) should also support that.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 9, 2023, 10:49pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/13 "2023-05-09T22:49:50Z")

</div>

> [@dlfivefifty](#):
>
> You can use `LowRankMatrix` in LowRankApprox.jl

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!

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 11, 2023, 11:58am UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/14 "2023-05-11T11:58:04Z")

</div>

> [@e3c6](#):
>
> Also, would it be a good idea to break `LowRankMatrix` apart into its own package?

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

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [May 11, 2023, 12:44pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/15 "2023-05-11T12:44:47Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 11, 2023, 12:53pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/16 "2023-05-11T12:53:58Z")

</div>

> [@jishnub](#):
>
> although what’s the point, instead of contributing changes to the package?

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.

---

<div class="post-metadata">

**Author:** ![jarl](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jarl/32/4365_2.png) [@jarl](https://discourse.julialang.org/u/jarl)\
**Post date:** [May 11, 2023, 2:10pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/17 "2023-05-11T14:10:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![e3c6](https://avatars.discourse-cdn.com/v4/letter/e/e79b87/32.png) [@e3c6](https://discourse.julialang.org/u/e3c6)\
**Post date:** [May 11, 2023, 3:03pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/18 "2023-05-11T15:03:40Z")

</div>

> [@jarl](#):
>
> since the data structure is not specific for low-rank matrices

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.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [July 7, 2023, 8:52pm UTC](https://discourse.julialang.org/t/low-rank-factorized-abstractmatrix/19838/19 "2023-07-07T20:52:02Z")

</div>

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