# ACAFact.jl

A simple and dependency-free adaptive cross approximation (ACA) factorization in Julia.

# Usage

Let's create a matrix with low numerical rank:

using LinearAlgebra
pts1 = range(0.0, 1.0, length=100)
pts2 = range(0.0, 1.0, length=110)
K = [exp(-abs2(pj - pk)) for pj in pts1, pk in pts2]
@show rank(K) # 10

One easy way to use `ACAFact.jl`

is to just provide your matrix and a maximum
allowed rank for the approximation:

for rk in (5, 10, 15, 20, 25)
(U, V) = aca(K, rk)
@show (rk, opnorm(K - U*V'))
end

Another thing you can do is provide a desired `opnorm`

error for `||K - U*V'||`

...but I wouldn't bet the farm that the error control is all that
guaranteed. Greedy deterministic pivoting schemes can be tricked, or at least
made very inefficient! It does work pretty well in toy applications at least
though, so it isn't useless.

for tol in (1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14)
(U, V) = aca(K, 100, tol=tol)
@show tol, size(U, 2), opnorm(K - U*V')
end

If you want something that is completely non-allocating, you can pre-allocate a
cache. Note that if you also provide a `tol`

here, you will terminate early but
the function won't remove the unused columns in `U`

and `V`

. The `rnk`

return
parameter gives you the rank of the approximation, so you should work with
`view(U, :, 1:rnk)`

, for example.

cache = ACAFact.ACACache(Float64, length(pts1), length(pts2), 50) # max rank 50
(U, V) = (zeros(length(pts1), 50), zeros(length(pts2), 50))
(rnk, _) = ACAFact.aca!(K, U, V, cache=cache) # zero allocations

One note here though: this `ACACache`

does carry state that gets used in the
factorizations, so if you want to factorize a new matrix that is the same size
as `K`

with the same cache, be sure to `ACAFact.resetcache!(cache)`

.

# Working with arbitrary matrix-like objects

To me, the best reason to use a greedy partially-pivoted factorization like the
ACA is to build low-rank approximations for matrices where you cannot afford to
pass over every entry even once. To accommodate this use case, `ACAFact.jl`

has
an interface

ACAFact.col!(buf, K, j)
ACAFact.row!(buf, K, j)

that expects `buf`

to be filled with the corresponding column/row of `K`

. So if
you have some cool operator that is defined implicitly or whatever, all you need
to do add special methods

ACAFact.col!(buf, K::MyCoolObject, j) = # ...
ACAFact.row!(buf, K::MyCoolObject, j) = # ...

and you can create low-rank approximations that never touch rows/columns of your
matrix that aren't selected as pivots. It would be nice if you could also
implement `Base.size(K::MyCoolObject)`

, but if for some reason you don't want to
you can compute the ACA of `K`

with `aca(K, [...], sz=(size_K_1, size_K_2))`

.