Say I have a very large vector or matrix M and I want to do M*M’ to obtain say a time by time matrix. This doesn’t work for huge matrices. I guess the issue is memory? Are there any common ways to solve this?
The answer depends on how you need it,
- if its to build an operator better build a SciMLOperator.jl matrix free operator (M M^T x) and solve what you need with Kyrlov
- if you need the matrix you will have to store it in disk and play one chunk at a time ie Mmap.jl.
- if you need some elements of the matrix you can go with LazyArrays.jl
- If you want something more “high-lever” i guess Dagger.jl + DiskArrays.jl / DiskArrayEngine.jl could do the trick but its very young and I cant comment on it.
I think you need to step back for a moment and think of the bigger picture.
First, if your matrix M fits in memory but your matrix MM^T does not, I’m guessing that M is a “tall” matrix (rows ≫ columns). Even though the textbook formula for your answer may involve MM^T, one typically does not explicitly form MM^T in practice — not only is the latter matrix huge, but it also squares the condition number (which may double the number of digits you lose to roundoff error). Typically one uses either QR or SVD factorizations of M instead.
(You can also form MM^T implicitly as a linear operator as @yolhan_mannes suggested above, then use an iterative solver, but this doesn’t solve the conditioning problem … and ill-conditioning may cause such an iterative solver to converge slowly.)
More generally, if even M itself is too big to fit in memory, then you usually want to look for structure in where the matrix comes from to let you implicitly construct it in terms of fast matrix–vector products. If M truly comes from unstructured data and doesn’t fit in memory, then you may be stuck with expensive “out-of-core” algorithms that involve reading M by chunks, but even then one often tries to shrink the data required by e.g. random sampling / “sketching”.
What problem are you trying to solve with MM^T? That is, what are you trying to do with this matrix? And how big is M and where does it come from?
Thanks for the thorough comments! Yeah, it’s a “tall” matrix as you mention. I want to try out different approaches that may be similar to what recurrence quantification analysis does (not necessarily the embedded delay addition part though), including approaches that probably should not be called like that because they don’t dichotomize the resulting matrices. And try different ways to construct the time by time matrix (self-, cross- similarity/distance). The problem is (lag-)synchrony estimation. M does fit in memory.
The point is that any algorithm that relies on constructing the whole Gram matrix G = MM^T will not scale, and you want to look for ways to skip over this to compute your final desired result implicitly without G. If you’re ultimately going to threshold G or something, for example, you might look into tree-based nearest-neighbor searches.
The very pragmatic answer often is: Then do try it out! (apologies for the trite generalities)
Create a bunch of smaller sample problems, try out your various algorithms, always plot whatever you want depending on N (including memory consumption and runtime and numerical errors).
You will see that M*M' scales quadratically with N for tall matrices, so this approach cannot be used for non-small data; but you can figure out all the other parts of your pipeline, and the issue is eminently solvable using various methods like @stevengj mentioned.
It is a typical mistake to do your experimentation with a single largish dataset: This means that one early problem (“gram matrix too big”) blocks your entire progress on the later parts that actually interest you.
On the other hand, if you already have later experimentations on a sequence of smaller problem sizes, then you know much better what you’re actually trying to do with M*M' and we can give much better advice on how to fix that.
(this approach doesn’t work for stuff with really bad scaling, where N+1 is a qualitatively different problem from N)
Thanks everyone for the excellent responses! I’m just trying to understand what would be some easy starting points (existing implementations etc.). For instance a basic one is taking the sum across diagonals in the gram matrix (which results in something akin to an autocorrelation curve, depending on implementation details).
If you simply need to access elements of the matrix, you can compute those entries as you need them. Or it’s very easy to make a custom type that does this (while pretending to be the full matrix) by using the array interface:
struct OuterProd{T, X<:AbstractVecOrMat{T}} <: AbstractMatrix{T}
x::X
end
Base.size(m::OuterProd) = (size(m.x, 1), size(m.x, 1))
Base.getindex(m::OuterProd, i::Int, j::Int) = m.x[j, :]' * m.x[i, :]
# if needed, could add methods for matrix multiplication, etc.
julia> x = ones(10^6, 10); m = OuterProd(x)
1000000×1000000 OuterProd{Float64, Matrix{Float64}}:
10.0 10.0 10.0 10.0 10.0 … 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 … 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
⋮ ⋱
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 … 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0
julia> sum(diag(m, 3)) # sum of 3rd superdiagonal
9.99997e6
Despite this matrix having 1 trillion elements, we compute them on-demand from X X^\top so it only requires the memory of X. You don’t pay for the cost of elements you don’t use.
If you’re going to multiply something by this matrix, you’ll want (in practice, need) to specialize multiplication to be efficient. By default it will access every element of this matrix (which is unecessarily expensive when it is large).
Maybe LazyArrays.jl has a version of this (with more features already baked in), but I couldn’t figure it out from the docs.
. That is insanely useful. Thanks for clarifying a part of the Julia documentation that was very obscure to me.
Not sure if I must create another thread now, but what if I want to get the sum across each diagonal? I tried this but there might be something better:
dind = [-(size(m,1)-1):0; 1:(size(m,1)-1)]
[sum(diag(m,x)) for x in dind]
It’s a hermitian matrix, so the +n diagonal is the conjugate of the -n diagonal (i.e., the same if real-valued). That can save you half the work.
Since a lot of the time here is spent in matrix multiplication, it’s probably more efficient to compute a column of X X^\top at a time and then accumulate those results into the proper diagonals (i.e., the first column contributes to the -n+1:0 diagonals, the second to the -n+2:+1 diagonals, etc). This can be combined with the “only compute half the diagonals” part above. This is more efficient because the underlying linear algebra functions can do one matrix-vector multiply more quickly than n vector-vector dot products. Unfortunately, the cute little wrapper I wrote needs another getindex method to do this efficiently (because I only implemented the single-element version).
Thanks. I’m asking with a XY’ in mind as well, so not just for symmetric matrices. Unfortunately I don’t understand the second paragraph. By the way, sum(diag(m, 3)) is very fast but things seem to get quite rather when I loop.
I have tried this, which probably is not the correct way to go since I’m rather new with this:
struct OuterProdDiag{T, X<:OuterProd{T}}
x::X
end
Base.size(m::OuterProdDiag) = size(m.x)
Base.getindex(m::OuterProdDiag, i::Int) = sum(diag(m.x,i))
o = OuterProdDiag(m)
I’m not sure whether more methods are needed here to get a curve with all the sums across diagonals.
Computing every diagonal is slow because it’s needing to compute every entry (of the massive matrix) to do it. While OuterProd solved the memory problem, there isn’t much to improve the computation part of it.
You don’t need a new type (struct) for this. I would just implement a function sumdiagonal(x, k) = sum(diagview(x, k)) (notice, it’s fine for any type of matrix).
And then something in the likes [sumdiagonal(m,k) for k in dind]? Not sure if I see the advantage over sum(diag(m,k)), at least for an 1000x1000 matrix.
julia> @btime [sumdiagonal(m,k) for k in dind];
22.197 ms (2002980 allocations: 91.61 MiB)
julia> @btime [sum(diag(m,k)) for k in dind];
22.886 ms (2008475 allocations: 100.00 MiB)
There’s no advantage. It’s doing the same thing. I suggested it because it looked like you were trying to do with OuterProdDiag. I didn’t see any reason for that to be a type so suggested it be a function instead.
But if you want it to be efficient, you need to make better use of the highly optimized linear algebra library. For example, using matrix-vector multiplies instead of (more) vector-vector dot products. Here’s a better version that does this and also only computes half the values because of symmetry:
# for comparison
sumdiagonals_simple(x) = [sum(diagview(x, d)) for d in -(size(x,1)-1):(size(x,2)-1)]
# more optimized
function sumdiagonals(A::OuterProd)
n = size(A, 1)
diags = zeros(eltype(A), 2n-1) # middle index (diagonal 0) is at index n
for i in 0:n-1
# compute X*X', one column (of the upper triangle) at a time
# then accumulate these values into diags
X = @view A.x[begin:begin+i, :] # the rows of the upper column
Xi = transpose(@view A.x[begin+i, :]) # one row of X (transpose so it's still a row)
diagv = @view diags[n+i:-1:n] # this upper column contributes to diagonals i:-1:0
mul!(diagv, X, Xi', true, true) # diagv .+= X * Xi'
end
diags[begin:n-1] = @view diags[end:-1:n+1] # reverse and copy to the other half
return diags
end
julia> m = OuterProd(ones(10^4, 3))
10000×10000 OuterProd{Float64, Matrix{Float64}}:
3.0 3.0 3.0 3.0 3.0 3.0 … 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 … 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
⋮ ⋮ ⋱ ⋮
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 … 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0
julia> sumdiagonals_simple(m) == sumdiagonals(m)
true
julia> @time sumdiagonals_simple(m);
2.270754 seconds (200.00 M allocations: 8.941 GiB, 14.84% gc time)
julia> @time sumdiagonals(m);
0.024738 seconds (6 allocations: 234.516 KiB)
I will leave the X Y^\top case to you for now (hint: it looks somewhat similar, though needs to compute full columns at a time).
Not sure what exactly you want to do, but just as a reference, there is this project computing RQA-based GitHub - EarthyScience/RQADeforestation.jl: Sentinel-1 based deforestation detection written in Julia · GitHub trend metrics. There the team also completely avoided allocating the whole recurrence matrix and did the whole linear regression on the fly by using online algorithms. In case you are interested @Fliks can add more details.
Interesting indeed!!
Thanks a million for this. I’m trying to understand it.
This was a mistake. I had the indexing backwards (it should have been diagv = @view diags[n+i:-1:n]). I edited my above post to fix it so make sure you grab the new version.
(That’s what I get for testing with X = ones(...). Should have used a random matrix.)
I got a lost with mul!(diagv, X, Xi', true, true). What are those trues about?