Fixed L-BFGS, included new code for BlockArray

src/AbstractOperators.jl

@@ -2,24 +2,24 @@ __precompile__()
 
 module AbstractOperators
 
-const RealOrComplex{T<:Real} = Union{T, Complex{T}}
-
 abstract type AbstractOperator end
 
 abstract type LinearOperator    <: AbstractOperator end
 abstract type NonLinearOperator <: AbstractOperator end
 
-import Base: A_mul_B!, Ac_mul_B! 
+import Base: A_mul_B!, Ac_mul_B!
 
 export LinearOperator,
        NonLinearOperator,
        AbstractOperator
 
-# deep stuff
+# Block stuff
+
+include("utilities/block.jl")
+#include("utilities/deep.jl") # TODO: remove this eventually
 
-include("utilities/deep.jl")
+# Predicates and properties
 
-# predicates and properties
 include("properties.jl")
 
 # Linear operators
@@ -43,7 +43,7 @@ include("linearoperators/Filt.jl")
 include("linearoperators/MIMOFilt.jl")
 include("linearoperators/Xcorr.jl")
 include("linearoperators/LBFGS.jl")
-include("linearoperators/BlkDiagLBFGS.jl")
+# include("linearoperators/BlkDiagLBFGS.jl")
 
 # Calculus rules
 
@@ -67,7 +67,4 @@ include("nonlinearoperators/SoftPlus.jl")
 # Syntax
 include("syntax.jl")
 
-
-
-
 end

src/linearoperators/LBFGS.jl

@@ -1,132 +1,125 @@
 export LBFGS, update!
 
-# TODO make Ac_mul_B!
-# Edit: Ac_mul_B! is not really needed for this operator
-# Edit2: you never known! anyway for completeness would be cool to have it!
 """
-`LBFGS(T::Type, dim::Tuple, Memory::Int)`
+`LBFGS(domainType::Type,dim_in::Tuple, M::Integer)`
 
-`LBFGS{N}(T::NTuple{N,Type}, dim::NTuple{N,Tuple}, M::Int)`
+`LBFGS(dim_in::Tuple, M::Integer)`
 
-`LBFGS(x::AbstractArray, Memory::Int)`
+`LBFGS(x::AbstractArray, M::Integer)`
 
 Construct a Limited-Memory BFGS `LinearOperator` with memory `M`. The memory of `LBFGS` can be updated using the function `update!`, where the current iteration variable and gradient (`x`, `grad`) and the previous ones (`x_prev` and `grad_prev`) are needed:
 
 ```
 julia> L = LBFGS(Float64,(4,),5)
 LBFGS  ℝ^4 -> ℝ^4
 
-julia> update!(L,x,x_prev,grad,grad_prev); #update memory
+julia> update!(L,x,x_prev,grad,grad_prev); # update memory
 
-julia> d = L*x;                            #compute new direction
+julia> d = L*grad; # compute new direction
 
 ```
-
 """
 
-mutable struct LBFGS{M, N, R <: Real, T <: Union{R, Complex{R}}, A<:AbstractArray{T,N}} <: LinearOperator
-	currmem::Int
-	curridx::Int
-	s::A
-	y::A
-	s_m::NTuple{M, A}
-	y_m::NTuple{M, A}
-	ys_m::Array{R, 1}
+mutable struct LBFGS{R, T <: BlockArray, M, I <: Integer} <: LinearOperator
+	currmem::I
+	curridx::I
+	s::T
+	y::T
+	s_M::Array{T, 1}
+	y_M::Array{T, 1}
+	ys_M::Array{R, 1}
 	alphas::Array{R, 1}
 	H::R
 end
 
-# Constructors
-#default
-function LBFGS(T::Type, dim::NTuple{N,Int}, M::Int) where {N}
-	s_m = tuple([deepzeros(T,dim) for i = 1:M]...)
-	y_m = tuple([deepzeros(T,dim) for i = 1:M]...)
-	s = deepzeros(T,dim)
-	y = deepzeros(T,dim)
-	R = real(T)
-	ys_m = zeros(R, M)
+#default constructor
+
+function LBFGS(domainType, dim_in, M::I) where {I <: Integer}
+	s_M = [blockzeros(domainType, dim_in) for i = 1:M]
+	y_M = [blockzeros(domainType, dim_in) for i = 1:M]
+	s = blockzeros(domainType, dim_in)
+	y = blockzeros(domainType, dim_in)
+	T = typeof(s)
+	R = typeof(domainType) <: Tuple  ? real(domainType[1]) : real(domainType)
+	ys_M = zeros(R, M)
 	alphas = zeros(R, M)
-	LBFGS{M,N,R,T,typeof(s)}(0, 0, s, y, s_m, y_m, ys_m, alphas, one(R))
+	LBFGS{R, T, M, I}(0, 0, s, y, s_M, y_M, ys_M, alphas, one(R))
 end
 
-LBFGS(x::AbstractArray,M::Int) = LBFGS(eltype(x),size(x),M)
+function LBFGS(dim_in, M::I) where {I <: Integer}
+	domainType = eltype(dim_in) <: Integer ? Float64 : ([Float64 for i in eachindex(dim_in)]...)
+	LBFGS(domainType, dim_in, M)
+end
+
+function LBFGS(x::T, M::I) where {T <: BlockArray, I <: Integer}
+	domainType = blockeltype(x)
+	dim_in = blocksize(x)
+	LBFGS(domainType, dim_in, M)
+end
 
 """
 `update!(L::LBFGS, x, x_prex, grad, grad_prev)`
 
-See `LBFGS` documentation.
-
+See the documentation for `LBFGS`.
 """
 
-function update!(L::LBFGS{M,N,R,T,A},
-		 x::A,
-		 x_prev::A,
-		 gradx::A,
-		 gradx_prev::A) where {M,N,R,T,A}
-
-	ys = update_s_y(L,x,x_prev,gradx,gradx_prev)
-
+function update!(L::LBFGS{R, T, M, I}, x::T, x_prev::T, gradx::T, gradx_prev::T) where {R, T, M, I}
+	L.s .= x .- x_prev
+	L.y .= gradx .- gradx_prev
+	ys = real(blockvecdot(L.s, L.y))
 	if ys > 0
 		L.curridx += 1
 		if L.curridx > M L.curridx = 1 end
 		L.currmem += 1
 		if L.currmem > M L.currmem = M end
-
-
-		yty = update_s_m_y_m(L,L.curridx)
-		L.ys_m[L.curridx] = ys
+		L.ys_M[L.curridx] = ys
+		blockcopy!(L.s_M[L.curridx], L.s)
+		blockcopy!(L.y_M[L.curridx], L.y)
+		yty = real(blockvecdot(L.y, L.y))
 		L.H = ys/yty
 	end
 	return L
 end
 
-function update_s_y(L::LBFGS{M,N,R,T,A}, x::A, x_prev::A, gradx::A, gradx_prev::A) where {M,N,R,T,A}
-	L.s .= (-).(x, x_prev)
-	L.y .= (-).(gradx, gradx_prev)
-	ys = real(vecdot(L.s,L.y))
-	return ys
-end
+# LBFGS operators are symmetric
 
-function update_s_m_y_m(L::LBFGS{M,N,R,T,A}, curridx::Int) where {M,N,R,T,A}
-	L.s_m[curridx] .=  L.s
-	L.y_m[curridx] .=  L.y
+Ac_mul_B!(x::T, L::LBFGS{R, T, M, I}, y::T) where {R, T, M, I} = A_mul_B!(x, L, y)
 
-	yty = real(vecdot(L.y,L.y))
-	return yty
-end
+# Two-loop recursion
 
-function A_mul_B!(d::A, L::LBFGS{M,N,R,T,A}, gradx::A) where {M,N,R,T,A}
-	d .= (-).(gradx)
+function A_mul_B!(d::T, L::LBFGS{R, T, M, I}, gradx::T) where {R, T, M, I}
+	d .= gradx
 	idx = loop1!(d,L)
 	d .= (*).(L.H, d)
 	d = loop2!(d,idx,L)
 end
 
-function loop1!(d::A, L::LBFGS{M,N,R,T,A}) where {M,N,R,T,A}
+function loop1!(d::T, L::LBFGS{R, T, M, I}) where {R, T, M, I}
 	idx = L.curridx
-	for i=1:L.currmem
-		L.alphas[idx] = real(vecdot(L.s_m[idx], d))/L.ys_m[idx]
-		d .-= L.alphas[idx].*L.y_m[idx]
+	for i = 1:L.currmem
+		L.alphas[idx] = real(blockvecdot(L.s_M[idx], d))/L.ys_M[idx]
+		d .-= L.alphas[idx] .* L.y_M[idx]
 		idx -= 1
 		if idx == 0 idx = M end
 	end
 	return idx
 end
 
-function loop2!(d::A, idx::Int, L::LBFGS{M,N,R,T,A}) where {M,N,R,T,A}
-	for i=1:L.currmem
+function loop2!(d::T, idx::Int, L::LBFGS{R, T, M, I}) where {R, T, M, I}
+	for i = 1:L.currmem
 		idx += 1
 		if idx > M idx = 1 end
-		beta = real(vecdot(L.y_m[idx], d))/L.ys_m[idx]
-		d .+= (L.alphas[idx].-beta).*L.s_m[idx]
+		beta = real(blockvecdot(L.y_M[idx], d))/L.ys_M[idx]
+		d .+= (L.alphas[idx] - beta) .* L.s_M[idx]
 	end
 	return d
 end
 
 # Properties
-  domainType(L::LBFGS{M,N,R,T,A}) where {M,N,R,T,A} = T
-codomainType(L::LBFGS{M,N,R,T,A}) where {M,N,R,T,A} = T
 
-size(A::LBFGS) = (size(A.s), size(A.s))
+domainType(L::LBFGS{R, T, M, I}) where {R, T, M, I} = blockeltype(L.y_M[1])
+codomainType(L::LBFGS{R, T, M, I}) where {R, T, M, I} = blockeltype(L.y_M[1])
+
+size(A::LBFGS) = (blocksize(A.s), blocksize(A.s))
 
 fun_name(A::LBFGS) = "LBFGS"