fixed blockvecdot for complex, added complex var test NonLinearCompose

Author: nantonel

src/utilities/block.jl

@@ -54,7 +54,8 @@ blockset!(y::Tuple, x) = blockset!.(y, x)
 blockset!(y::AbstractArray, x) = (y .= x)
 
 blockvecdot(x::T, y::T) where {T <: Tuple} = sum(blockvecdot.(x,y))
-blockvecdot(x::AbstractArray{R}, y::AbstractArray{R}) where {R <: Number} = vecdot(x, y)
+blockvecdot(x::AbstractArray{R}, y::AbstractArray{R}) where {R <: Number} = real(vecdot(x, y))
+# inner product must be always real see section 4.2 of TFOCS manual
 
 blockzeros(t::Tuple, s::Tuple) = blockzeros.(t, s)
 blockzeros(t::Type, n::NTuple{N, Integer} where {N}) = zeros(t, n)

test/test_block.jl

@@ -48,8 +48,8 @@ blockset!(yb,xb)
 
 z = blockvecdot(x,x2)
 zb = blockvecdot(xb,x2b)
-@test z == vecdot(x,x2)
-@test zb == vecdot(xb[1],x2b[1])+vecdot(xb[2],x2b[2])+vecdot(xb[3],x2b[3])
+@test z == real(vecdot(x,x2))
+@test zb == real(vecdot(xb[1],x2b[1])+vecdot(xb[2],x2b[2])+vecdot(xb[3],x2b[3]))
 
 y = blockzeros(x)
 yb = blockzeros(xb)

test/test_nonlinear_operators_calculus.jl

@@ -180,7 +180,7 @@ y, grad = test_NLop(op,x,r,verb)
 Y = A*x+opB*x
 @test vecnorm(Y-y) <1e-8
 
-#testing NonLinearCompose
+## testing NonLinearCompose
 
 #with vectors
 n,m = 3,4
@@ -208,7 +208,25 @@ Y = A*x[1]*B*x[2]
 @test norm(Y - y) <= 1e-12
 
 #further test on gradient with analytical formulas
-grad2 =  ((B*x[2])*(A'*r)')', B'*(A*x[1])'*r
+grad2 =  (A'*r)*(B*x[2])', B'*(A*x[1])'*r
+@test vecnorm(grad[1]-grad2[1]) <1e-7
+@test vecnorm(grad[2]-grad2[2]) <1e-7
+
+#with complex matrices
+l,m1,m2,n1,n2 = 2,3,4,5,6
+x = (randn(m1,m2)+im*randn(m1,m2),randn(n1,n2)+im*randn(n1,n2))
+A = randn(l,m1) +im*randn(l,m1)
+B = randn(m2,n1)+im*randn(m2,n1)
+r = randn(l,n2) +im*randn(l,n2)
+
+P = NonLinearCompose( MatrixOp(A,m2), MatrixOp(B,n2) )
+y, grad = test_NLop(P,x,r,verb)
+
+Y = A*x[1]*B*x[2]
+@test norm(Y - y) <= 1e-12
+
+##test on gradient with analytical formulas
+grad2 =  (A'*r)*(B*x[2])', B'*(A*x[1])'*r
 @test vecnorm(grad[1]-grad2[1]) <1e-7
 @test vecnorm(grad[2]-grad2[2]) <1e-7
 

test/utils.jl

@@ -26,7 +26,7 @@ function test_op(A::AbstractOperator, x, y, verb::Bool = false)
   return Ax
 end
 
-########### Test for LinearOperators
+########### Test for NonLinearOperators
 function test_NLop(A::AbstractOperator, x, y, verb::Bool = false)
 
 	verb && (println(),println(A))
@@ -54,9 +54,11 @@ function test_NLop(A::AbstractOperator, x, y, verb::Bool = false)
 
 	@test AbstractOperators.blockvecnorm(grad .- grad2) < 1e-8
 
-	grad3 = gradient_fd(A,Ax,x,y) #calculate gradient using finite differences
+	if all(isreal.(grad))  # currently finite difference gradient not working with complex variables 
+		grad3 = gradient_fd(A,Ax,x,y) #calculate gradient using finite differences
 
-	@test AbstractOperators.blockvecnorm(grad .- grad3) < 1e-4
+		@test AbstractOperators.blockvecnorm(grad .- grad3) < 1e-4
+	end
 
 	return Ax, grad
 end