With current TensorProductOperator implementation,

Lines 684 to 724 in c948888

    
           function LinearAlgebra.ldiv!(v::AbstractVecOrMat, L::TensorProductOperator, u::AbstractVecOrMat) 
        
               @assert L.isset "cache needs to be set up for operator of type $(typeof(L)). 
        
               set up cache by calling cache_operator(L::AbstractSciMLOperator, u::AbstractArray)" 
        
               mi, ni = size(L.inner) 
        
               mo, no = size(L.outer) 
        
               k = size(u, 2) 
        
               perm = (2, 1, 3) 
        
               C1, C2, C3, _ = L.cache 
        
               U = _reshape(u, (ni, no*k)) 
        
               """ 
        
                   v .= kron(B, A) ldiv u 
        
                   V .= (A ldiv U) / B' 
        
               """ 
        
               # C .= A \ U 
        
               ldiv!(C1, L.inner, U) 
        
               # V .= C / B' <===> V' .= B \ C' 
        
               if k>1 
        
                   if L.outer isa IdentityOperator 
        
                       copyto!(v, C1) 
        
                   else 
        
                       C1 = _reshape(C1, (mi, no, k)) 
        
                       permutedims!(C2, C1, perm) 
        
                       C2 = _reshape(C2, (no, mi*k)) 
        
                       ldiv!(C3, L.outer, C2) 
        
                       C3 = _reshape(C3, (mo, mi, k)) 
        
                       V  = _reshape(v , (mi, mo, k)) 
        
                       permutedims!(V, C3, perm) 
        
                   end 
        
               else 
        
                   V  = _reshape(v , (mi, mo)) 
        
                   C1 = _reshape(C1, (mi, no)) 
        
                   ldiv!(transpose(V), L.outer, transpose(C1)) 
        
               end 
        
               v 
        
           end

reference computation: - a single component matvec.

using LinearAlgebra, BenchmarkTools

u = rand(12, 12 * 100)
v = rand(12, 12 * 100)

A = rand(12, 12)

@btime mul!($v, $A, $u); # 23.833 μs (0 allocations: 0 bytes)

Id = SciMLOperators.IdentityOperator{12}()
@btime mul!($v, $Id, $u); # 4.822 μs (0 allocations: 0 bytes)

We want to be as close to 24+5 = 30\mu s for both I \otimes D and D \otimes I. Once we make TensorProducts efficient, we can figure out how to write a separate TensorSum wrapper later.

Other component operations:

permutedims!

using LinearAlgebra, BenchmarkTools
u = rand(12, 12, 100)
v = copy(u)
perm = (2, 1, 3)
@btime permutedims!($v, $u, $perm) # 8.980 μs (0 allocations: 0 bytes)

mul! with identity - doing good over here

using SciMLOperators, LinearAlgebra, BenchmarkTools
u = rand(12, 12* 100)
v = copy(u)
Id = SciMLOperators.IdentityOperator{12}()
@btime mul!($v, $Id, $u); #  4.790 μs (0 allocations: 0 bytes)
@btime mul!($v, I, $u); # 11.729 μs (0 allocations: 0 bytes)

I \otimes D - this is pretty good. nice and close to reference computation.

using SciMLOperators, LinearAlgebra, BenchmarkTools

u = rand(12^2, 100)
v = rand(12^2, 100)

A = SciMLOperators.IdentityOperator{12}() # outer
B = rand(12, 12) # inner

T = TensorProductOperator(A, B)
T = cache_operator(T, u)

@btime mul!($v, $T, $u); # 28.846 μs (1 allocation: 32 bytes)

D \otimes I - this is not so hot. Let's also try to this vector-wise

using SciMLOperators, LinearAlgebra, BenchmarkTools

u = rand(12^2, 100)
v = rand(12^2, 100)

A = rand(12, 12) # outer
B = SciMLOperators.IdentityOperator{12}() # inner

T = TensorProductOperator(A, B)
T = cache_operator(T, u)

@btime mul!($v, $T, $u); # 63.360 μs (7 allocations: 304 bytes)

# this is not too hot, let's try doing it vector-wise. 
u = [rand(12*12) for i=1:100]
v = [rand(12*12) for i=1:100]

T = cache_operator(T, u[1])

@btime for i=1:length($u) mul!($v[i], $T, $u[i]) end; # 228.426 μs (200 allocations: 2.00 MiB)

# this is worse :/

Batchwise component matvec

julia> using SciMLOperators, LinearAlgebra, BenchmarkTools; A = rand(12,12); u = [rand(12) for i=1:12*100]; v = copy(u); @btime for i=1:length($u); mul!($v[i], $A, $u[i]); end;
  104.498 μs (0 allocations: 0 bytes) # let's just not do batchwise

	function LinearAlgebra.ldiv!(v::AbstractVecOrMat, L::TensorProductOperator, u::AbstractVecOrMat)
	@assert L.isset "cache needs to be set up for operator of type $(typeof(L)).
	set up cache by calling cache_operator(L::AbstractSciMLOperator, u::AbstractArray)"

	mi, ni = size(L.inner)
	mo, no = size(L.outer)
	k = size(u, 2)

	perm = (2, 1, 3)
	C1, C2, C3, _ = L.cache
	U = _reshape(u, (ni, no*k))

	"""
	v .= kron(B, A) ldiv u
	V .= (A ldiv U) / B'
	"""

	# C .= A \ U
	ldiv!(C1, L.inner, U)

	# V .= C / B' <===> V' .= B \ C'
	if k>1
	if L.outer isa IdentityOperator
	copyto!(v, C1)
	else
	C1 = _reshape(C1, (mi, no, k))
	permutedims!(C2, C1, perm)
	C2 = _reshape(C2, (no, mi*k))
	ldiv!(C3, L.outer, C2)
	C3 = _reshape(C3, (mo, mi, k))
	V = _reshape(v , (mi, mo, k))
	permutedims!(V, C3, perm)
	end
	else
	V = _reshape(v , (mi, mo))
	C1 = _reshape(C1, (mi, no))
	ldiv!(transpose(V), L.outer, transpose(C1))
	end

	v
	end

Uh oh!

Uh oh!

[later] kronsum #53

Description

Metadata

Metadata

Assignees

Labels

Type

Fields

Projects

Milestone

Relationships

Development

Issue actions