Merge pull request #12 from nantonel/master

Fix to nonlinear displaced operators

Author: nantonel

benchmarks/LineSpectraEstimation.jl

@@ -121,12 +121,12 @@ end
 
 #StructuredOptimization non-Matrix Free
 function solve_problem!(slv::S, x0, y, K, F::A, Fc, lambda, lambda_m) where {S <: StructuredOptimization.ForwardBackwardSolver, A <: AbstractMatrix}
-	it, = @minimize ls(F*x0-y)+lambda_m*norm(x0,1) with slv
+    it, = @minimize ls(F*x0-complex(y))+lambda_m*norm(x0,1) with slv
 	return x0, it
 end
 
 function solve_problem_ncvx!(slv::S, x0, y, K, F::A, Fc, lambda, lambda_m) where {S <: StructuredOptimization.ForwardBackwardSolver, A <: AbstractMatrix}
-	it, = @minimize ls(F*x0-y) st norm(x0,0) <= 2*K with slv
+    it, = @minimize ls(F*x0-complex(y)) st norm(x0,0) <= 2*K with slv
 	return x0, it
 end
 

benchmarks/run_demos.jl

@@ -15,9 +15,9 @@
 #results = MatrixDecomposition.run_demo()
 #MatrixDecomposition.show_results(results...)
 
-#include("DNN.jl")
-#results = DNN.run_demo()
-#DNN.show_results(results...)
+include("DNN.jl")
+results = DNN.run_demo()
+DNN.show_results(results...)
 
 #include("TotalVariation.jl")
 #results = TotalVariation.run_demo()

docs/src/expressions.md

@@ -75,6 +75,12 @@ variation
 
 ### Nonlinear mappings
 ```@docs
+sin
+cos
+atan
+tanh
+exp
+pow
 sigmoid
 ```
 
@@ -86,5 +92,6 @@ Notice that these commands work also for the `Term`s described in [Functions and
 ```@docs
 variables
 operator
-displacement
+affine
+AbstractOperators.displacement
 ```

docs/src/tutorial.md

@@ -2,7 +2,7 @@
 
 ## Standard problem formulation
 
-Currently with `StructuredOptimization.jl` you can solve problems of the form
+Currently with StructuredOptimization.jl one can solve problems of the form
 
 ```math
 \underset{ \mathbf{x} }{\text{minimize}} \ f(\mathbf{x}) + g(\mathbf{x}),
@@ -18,7 +18,7 @@ The *least absolute shrinkage and selection operator* (LASSO) belongs to this cl
 \underset{ \mathbf{x} }{\text{minimize}} \ \tfrac{1}{2} \| \mathbf{A} \mathbf{x} - \mathbf{y} \|^2+ \lambda \| \mathbf{x} \|_1.
 ```
 
-Here the squared norm $\tfrac{1}{2} \| \mathbf{A} \mathbf{x} - \mathbf{y} \|^2$ is a *smooth* function $f$ wherelse the $l_1$-norm is a *nonsmooth* function $g$. This problem can be solved with only few lines of code:
+Here the squared norm $\tfrac{1}{2} \| \mathbf{A} \mathbf{x} - \mathbf{y} \|^2$ is a *smooth* function $f$ whereas the $l_1$-norm is a *nonsmooth* function $g$. This problem can be solved with only few lines of code:
 
 ```julia
 julia> using StructuredOptimization
@@ -52,7 +52,7 @@ julia> ~x                             # inspect solution
 
 It is possible to access to the solution by typing `~x`.
 By default variables are initialized by `Array`s of zeros.
-Different initializations can be set during construction `x = Variable( [1.; 0.; ...] )` or by assignement `~x .= [1.; 0.; ...]`.
+Different initializations can be set during construction `x = Variable( [1.; 0.; ...] )` or by assignment `~x .= [1.; 0.; ...]`.
 
 ## Constrained optimization
 
@@ -70,12 +70,12 @@ can be converted into an *indicator function*
 ```math
 g(\mathbf{x}) = \delta_{\mathcal{S}} (\mathbf{x}) =  \begin{cases}
     0       & \text{if} \ \mathbf{x} \in \mathcal{S},\\
-    +\infty & \text{otherwise},
+    +\infty & \text{otherwise}.
     \end{cases}
 ```
 
-to obtain the standard form. Constraints are treated as *nonsmooth functions*.
-This conversion is automatically performed by `StructuredOptimization.jl`.
+Constraints are treated as *nonsmooth functions*.
+This conversion is automatically performed by StructuredOptimization.jl.
 For example, the non-negative deconvolution problem:
 
 ```math
@@ -85,7 +85,7 @@ For example, the non-negative deconvolution problem:
 \end{align*}
 ```
 
-where $*$ stands fof convoluton and $\mathbf{h}$ contains the taps of a finite impluse response,
+where $*$ stands for convolution and $\mathbf{h}$ contains the taps of a finite impulse response filter,
 can be solved using the following lines of code:
 
 ```julia
@@ -102,7 +102,7 @@ julia> @minimize ls(conv(x,h)-y) st x >= 0.
 !!! note
 
     The convolution mapping was applied to the variable `x` using `conv`.
-    `StructuredOptimization.jl` provides a set of functions that can be
+    StructuredOptimization.jl provides a set of functions that can be
     used to apply specific operators to variables and create mathematical
     expression. The available functions can be found in [Mappings](@ref).
     In general it is more convenient to use these functions instead of matrices,
@@ -134,7 +134,7 @@ julia> @minimize ls(X1*X2-Y) st X1 >= 0., X2 >= 0.
 
 ## Limitations
 
-Currently `StructuredOptimization.jl` supports only *proximal gradient algorithms* (i.e., *forward-backward splitting* base), which require specific properties of the nonsmooth functions and costraint to be applicable. In particular, the nonsmooth functions must have an *efficiently computable proximal mapping*.
+Currently StructuredOptimization.jl supports only *proximal gradient algorithms* (i.e., *forward-backward splitting* base), which require specific properties of the nonsmooth functions and constraint to be applicable. In particular, the nonsmooth functions must have an *efficiently computable proximal mapping*.
 
 If we express the nonsmooth function $g$ as the composition of
 a function $\tilde{g}$ with a linear operator $A$:
@@ -148,7 +148,7 @@ then the proximal mapping of $g$ is efficiently computable if either of the foll
 
 1. Operator $A$ is a *tight frame*, namely it satisfies $A A^* = \mu Id$, where $\mu \geq 0$, $A^*$ is the adjoint of $A$, and $Id$ is the identity operator.
 
-2. Function $g$ is the *separable sum* $g(\mathbf{x}) = \sum_j h_j (B_j \mathbf{x}_j)$, where $\mathbf{x}_j$ are non-overlapping slices of $\mathbf{x}$, and $B_j$ are tight frames.
+2. Function $g$ is a *separable sum* $g(\mathbf{x}) = \sum_j h_j (B_j \mathbf{x}_j)$, where $\mathbf{x}_j$ are non-overlapping slices of $\mathbf{x}$, and $B_j$ are tight frames.
 
 Let us analyze these rules with a series of examples.
 The LASSO example above satisfy the first rule:
@@ -157,8 +157,8 @@ The LASSO example above satisfy the first rule:
 julia> @minimize ls( A*x - y ) + λ*norm(x, 1)
 ```
 
-since the non-smooth function $\lambda \| \cdot \|_1$ is not composed with any operator (or equivalently is composed with $Id$ which is a tight frame).
-Also the following problem would be accepted:
+since the nonsmooth function $\lambda \| \cdot \|_1$ is not composed with any operator (or equivalently is composed with $Id$ which is a tight frame).
+Also the following problem would be accepted by StructuredOptimization.jl:
 
 ```julia
 julia> @minimize ls( A*x - y ) + λ*norm(dct(x), 1)

src/solvers/build_solve.jl

@@ -39,11 +39,13 @@ function build(terms::Tuple, solver::ForwardBackwardSolver)
 			append!(kwargs, [(:Aq, Aq)])
 		end
 		if !isempty(smooth)
-			fs = extract_functions(smooth)
-			As = extract_operators(x, smooth)
 			if is_linear(smooth)
+                fs = extract_functions(smooth)
+                As = extract_operators(x, smooth)
 				append!(kwargs, [(:As, As)])
 			else
+                fs = extract_functions_nodisp(smooth)
+                As = extract_affines(x, smooth)
 				fs = PrecomposeNonlinear(fs, As)
 			end
 			append!(kwargs, [(:fs, fs)])

src/solvers/terms_extract.jl

@@ -24,6 +24,14 @@ end
 extract_functions{N}(t::NTuple{N,Term}) = SeparableSum(extract_functions.(t))
 extract_functions(t::Tuple{Term}) = extract_functions(t[1])
 
+# extract functions from terms without displacement
+function extract_functions_nodisp(t::Term)
+	f = t.lambda == 1. ? t.f : Postcompose(t.f, t.lambda)
+	return f
+end
+extract_functions_nodisp{N}(t::NTuple{N,Term}) = SeparableSum(extract_functions_nodisp.(t))
+extract_functions_nodisp(t::Tuple{Term}) = extract_functions_nodisp(t[1])
+
 # extract operators from terms
 
 # returns all operators with an order dictated by xAll
@@ -43,6 +51,48 @@ function extract_operators{N,M}(xAll::NTuple{N,Variable}, t::NTuple{M,Term})
 	return vcat(ops...)
 end
 
+sort_and_extract_operators(xAll::Tuple{Variable}, t::Term) = operator(t)
+
+function sort_and_extract_operators{N}(xAll::NTuple{N,Variable}, t::Term)
+	p = zeros(Int,N)
+	xL = variables(t)
+	for i in eachindex(xAll)
+		p[i] = findfirst( xi -> xi == xAll[i], xL)
+	end
+	return operator(t)[p]
+end
+
+# extract affines from terms
+
+# returns all affines with an order dictated by xAll
+
+#single term, single variable
+extract_affines(xAll::Tuple{Variable}, t::Term)  = affine(t)
+
+extract_affines{N}(xAll::NTuple{N,Variable}, t::Term) = extract_affines(xAll, (t,))
+
+#multiple terms, multiple variables
+function extract_affines{N,M}(xAll::NTuple{N,Variable}, t::NTuple{M,Term})
+	ops = ()
+	for ti in t
+		tex = expand(xAll,ti)
+		ops = (ops...,sort_and_extract_affines(xAll,tex))
+	end
+	return vcat(ops...)
+end
+
+sort_and_extract_affines(xAll::Tuple{Variable}, t::Term) = affine(t)
+
+function sort_and_extract_affines{N}(xAll::NTuple{N,Variable}, t::Term)
+	p = zeros(Int,N)
+	xL = variables(t)
+	for i in eachindex(xAll)
+		p[i] = findfirst( xi -> xi == xAll[i], xL)
+	end
+	return affine(t)[p]
+end
+
+# expand term domain dimensions
 function expand{N,T1,T2,T3}(xAll::NTuple{N,Variable}, t::Term{T1,T2,T3})
 	xt   = variables(t)
 	C    = codomainType(operator(t))
@@ -57,17 +107,6 @@ function expand{N,T1,T2,T3}(xAll::NTuple{N,Variable}, t::Term{T1,T2,T3})
 	return Term(t.lambda, t.f, ex)
 end
 
-sort_and_extract_operators(xAll::Tuple{Variable}, t::Term) = operator(t)
-
-function sort_and_extract_operators{N}(xAll::NTuple{N,Variable}, t::Term)
-	p = zeros(Int,N)
-	xL = variables(t)
-	for i in eachindex(xAll)
-		p[i] = findfirst( xi -> xi == xAll[i], xL)
-	end
-	return operator(t)[p]
-end
-
 # extract function and merge operator
 function extract_merge_functions(t::Term)
     if is_sliced(t)

src/syntax/expressions/abstractOperator_bind.jl

@@ -20,12 +20,7 @@ julia> reshape(A*x-b,2,5)
 function reshape(a::AbstractExpression, dims...)
 	A = convert(Expression,a)
 	op = Reshape(A.L, dims...)
-	if typeof(displacement(A)) <: Number
-		d = displacement(A)
-	else
-		d = reshape(displacement(A), dims...)
-	end
-	return Expression{length(A.x)}(A.x,op,d)
+	return Expression{length(A.x)}(A.x,op)
 end
 #Reshape
 
@@ -39,6 +34,11 @@ imported = [:getindex :GetIndex;
             :conv     :Conv;
             :xcorr    :Xcorr;
             :filt     :Filt;
+            :exp      :Exp;
+            :cos      :Cos;
+            :sin      :Sin;
+            :atan     :Atan;
+            :tanh     :Tanh;
            ]
 
 exported = [:finitediff :FiniteDiff;
@@ -47,6 +47,7 @@ exported = [:finitediff :FiniteDiff;
             :zeropad    :ZeroPad;
             :sigmoid    :Sigmoid;
             :σ          :Sigmoid; #alias
+            :pow        :Pow; #alias
            ]
 
 #importing functions from Base
@@ -396,3 +397,75 @@ See documentation of `AbstractOperator.Sigmoid`.
 """
 sigmoid    
 σ          
+
+"""
+`exp(x::AbstractExpression)`
+
+Exponential function:
+```math
+e^{ \\mathbf{x}  }
+```
+
+See documentation of `AbstractOperator.Exp`.  
+"""
+exp 
+
+"""
+`sin(x::AbstractExpression)`
+
+Sine function:
+```math
+\\sin( \\mathbf{x}  )
+```
+
+See documentation of `AbstractOperator.Sin`.  
+"""
+sin 
+
+"""
+`cos(x::AbstractExpression)`
+
+Cosine function:
+```math
+\\cos( \\mathbf{x}  )
+```
+
+See documentation of `AbstractOperator.Cos`.  
+"""
+cos 
+
+"""
+`atan(x::AbstractExpression)`
+
+Inverse tangent function:
+```math
+\\tan^{-1}( \\mathbf{x}  )
+```
+
+See documentation of `AbstractOperator.Atan`.  
+"""
+atan 
+
+"""
+`tanh(x::AbstractExpression)`
+
+Hyperbolic tangent function:
+```math
+\\tanh ( \\mathbf{x}  )
+```
+
+See documentation of `AbstractOperator.Tanh`.  
+"""
+tanh 
+
+"""
+`pow(x::AbstractExpression, n)`
+
+Elementwise power 'n' of 'x':
+```math
+x_i^{n} \\ \\forall \\  i = 0,1, \\dots
+```
+
+See documentation of `AbstractOperator.Pow`.  
+"""
+pow