diff --git a/demos/critical_filtering.py b/demos/critical_filtering.py
index de8108f4d77dd4f462733ed37c48816c30d8e03c..e5f88adfea5b03bc345792f28c6bfb08d1adfedf 100644
--- a/demos/critical_filtering.py
+++ b/demos/critical_filtering.py
@@ -34,7 +34,7 @@ if __name__ == "__main__":
     # Instrument._diagonal.val[64:512-64, 64:512-64] = 0
 
     # Add a harmonic transformation to the instrument
-    R = ift.ComposedOperator([fft, Instrument])
+    R = Instrument*fft
 
     noise = 1.
     N = ift.DiagonalOperator(ift.Field.full(s_space, noise).weight(1))
diff --git a/demos/log_normal_wiener_filter.py b/demos/log_normal_wiener_filter.py
index b6bd4f6046fd4ff755523892a7d7f2bdc8e93cff..7e54d6ba10f75f453309e79429c30555d40e4d1f 100644
--- a/demos/log_normal_wiener_filter.py
+++ b/demos/log_normal_wiener_filter.py
@@ -36,7 +36,7 @@ if __name__ == "__main__":
     R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
                              exposure=(mask,))
     data_domain = R.target[0]
-    R_harmonic = ift.ComposedOperator([fft, R])
+    R_harmonic = R*fft
 
     # Setting up the noise covariance and drawing a random noise realization
     ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
@@ -66,20 +66,16 @@ if __name__ == "__main__":
     # m3 = fft(me3[0].position)
 
     # Plotting
-    ift.plotting.plot(mock_signal, name='mock_signal.png', colormap="plasma",
-                      xlabel="Pixel Index", ylabel="Pixel Index")
+    plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
+                "colormap": "Planck-like"}
+    ift.plotting.plot(mock_signal, name="mock_signal.png", **plotdict)
     logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
     ift.plotting.plot(ift.Field(signal_space,
                                 val=ift.dobj.from_global_data(logdata)),
-                      name="log_of_data.png", colormap="plasma",
-                      xlabel="Pixel Index", ylabel="Pixel Index")
-    # ift.plotting.plot(m1,name='m_LBFGS.png', colormap="plasma",
-    #                   xlabel="Pixel Index", ylabel="Pixel Index")
-    ift.plotting.plot(m2, name='m_Newton.png', colormap="plasma",
-                      xlabel="Pixel Index", ylabel="Pixel Index")
-    # ift.plotting.plot(m3, name='m_SteepestDescent.png',
-    #                   colormap="plasma", xlabel="Pixel Index",
-    #                   ylabel="Pixel Index")
+                      name="log_of_data.png", **plotdict)
+    # ift.plotting.plot(m1,name='m_LBFGS.png', **plotdict)
+    ift.plotting.plot(m2, name='m_Newton.png', **plotdict)
+    # ift.plotting.plot(m3, name='m_SteepestDescent.png', **plotdict)
 
     # Probing the variance
     class Proby(ift.DiagonalProberMixin, ift.Prober):
@@ -89,4 +85,4 @@ if __name__ == "__main__":
 
     sm = ift.FFTSmoothingOperator(signal_space, sigma=0.02)
     variance = sm(proby.diagonal.weight(-1))
-    ift.plotting.plot(variance, name='variance.png')
+    ift.plotting.plot(variance, name='variance.png', **plotdict)
diff --git a/demos/nonlinear_critical_filter.py b/demos/nonlinear_critical_filter.py
index ae0059575a2a94a8b814f61ea499d439dd81a84b..0f7ecc0756d2cd8c7ccff94f76aaf4d4de7657be 100644
--- a/demos/nonlinear_critical_filter.py
+++ b/demos/nonlinear_critical_filter.py
@@ -50,8 +50,8 @@ if __name__ == "__main__":
 
     MaskOperator = ift.DiagonalOperator(mask)
     InstrumentResponse = ift.ResponseOperator(s_space, sigma=[0.0])
-    MeasurementOperator = ift.ComposedOperator([MaskOperator,
-                                                InstrumentResponse])
+    MeasurementOperator = InstrumentResponse*MaskOperator
+
     d_space = MeasurementOperator.target
 
     noise_covariance = ift.Field(d_space, val=noise_level**2).weight()
diff --git a/demos/paper_demos/cartesian_wiener_filter.py b/demos/paper_demos/cartesian_wiener_filter.py
index 655d7cf015caab288107dbb4427496064be342b9..95468b64aa9fe8b00e9717aabd6cdf90dd894045 100644
--- a/demos/paper_demos/cartesian_wiener_filter.py
+++ b/demos/paper_demos/cartesian_wiener_filter.py
@@ -44,7 +44,7 @@ if __name__ == "__main__":
 
     mock_power_2 = ift.PS_field(power_space_2, power_spectrum_2)
 
-    fft = ift.ComposedOperator((fft_1, fft_2))
+    fft = fft_2*fft_1
 
     mock_power = ift.Field(
         (power_space_1, power_space_2),
@@ -74,7 +74,7 @@ if __name__ == "__main__":
                              sigma=(response_sigma_1, response_sigma_2),
                              exposure=(mask_1, mask_2))
     data_domain = R.target
-    R_harmonic = ift.ComposedOperator([fft, R])
+    R_harmonic = R*fft
 
     # Setting up the noise covariance and drawing a random noise realization
     ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
@@ -105,13 +105,15 @@ if __name__ == "__main__":
 
     plot_space = ift.RGSpace((N_pixels_1, N_pixels_2))
     sm = ift.FFTSmoothingOperator(plot_space, sigma=0.03)
+    plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
+                "colormap": "Planck-like"}
     ift.plotting.plot(
         ift.log(ift.sqrt(sm(ift.Field(plot_space, val=variance.val.real)))),
         name='uncertainty.png', zmin=0., zmax=3., title="Uncertainty map",
-        colormap="Planck-like")
+        **plotdict)
     ift.plotting.plot(ift.Field(plot_space, val=mock_signal.val.real),
-                      name='mock_signal.png', colormap="Planck-like")
+                      name='mock_signal.png', **plotdict)
     ift.plotting.plot(ift.Field(plot_space, val=data.val.real),
-                      name='data.png', colormap="Planck-like")
+                      name='data.png', **plotdict)
     ift.plotting.plot(ift.Field(plot_space, val=m.val.real),
-                      name='map.png', colormap="Planck-like")
+                      name='map.png', **plotdict)
diff --git a/demos/paper_demos/wiener_filter.py b/demos/paper_demos/wiener_filter.py
index 716d469107b424a632a0ec808760462bcff0c556..9e71c50ba6975e0d1e68d37102f90c2ade31b3d3 100644
--- a/demos/paper_demos/wiener_filter.py
+++ b/demos/paper_demos/wiener_filter.py
@@ -38,7 +38,7 @@ if __name__ == "__main__":
     R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
                              exposure=(mask,))
     data_domain = R.target[0]
-    R_harmonic = ift.ComposedOperator([fft, R])
+    R_harmonic = R * fft
 
     # Setting up the noise covariance and drawing a random noise realization
     ndiag = ift.Field.full(data_domain, mock_signal.var()/signal_to_noise)
@@ -67,12 +67,10 @@ if __name__ == "__main__":
     variance = ift.sqrt(sm(proby.diagonal.weight(-1)))
 
     # Plotting
-    ift.plotting.plot(variance, name="uncertainty.png", xlabel='Pixel index',
-                      ylabel='Pixel index')
-    ift.plotting.plot(mock_signal, name="mock_signal.png",
-                      xlabel='Pixel index', ylabel='Pixel index')
+    plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
+                "colormap": "Planck-like"}
+    ift.plotting.plot(variance, name="uncertainty.png", **plotdict)
+    ift.plotting.plot(mock_signal, name="mock_signal.png", **plotdict)
     ift.plotting.plot(ift.Field(signal_space, val=data.val),
-                      name="data.png", xlabel='Pixel index',
-                      ylabel='Pixel index')
-    ift.plotting.plot(m, name="map.png", xlabel='Pixel index',
-                      ylabel='Pixel index')
+                      name="data.png", **plotdict)
+    ift.plotting.plot(m, name="map.png", **plotdict)
diff --git a/demos/wiener_filter_via_curvature.py b/demos/wiener_filter_via_curvature.py
index d79c837f910a1b375473613b34849b2af1c46958..57383f342b183d317c313032c2baf401f9825b8d 100644
--- a/demos/wiener_filter_via_curvature.py
+++ b/demos/wiener_filter_via_curvature.py
@@ -53,7 +53,7 @@ if __name__ == "__main__":
     R = ift.ResponseOperator(signal_space, sigma=(response_sigma,),
                              exposure=(exposure,))
     data_domain = R.target[0]
-    R_harmonic = ift.ComposedOperator([fft, R])
+    R_harmonic = R*fft
 
     N = ift.DiagonalOperator(
         ift.Field.full(data_domain,
diff --git a/demos/wiener_filter_via_hamiltonian.py b/demos/wiener_filter_via_hamiltonian.py
index 5676d2ef330b9530a89bd9bf7e0be8217f9ecf11..1a0b205f3b6b289f7b10c553e18e49e53cfadcfb 100644
--- a/demos/wiener_filter_via_hamiltonian.py
+++ b/demos/wiener_filter_via_hamiltonian.py
@@ -18,7 +18,7 @@ if __name__ == "__main__":
 
     # Choosing the prior correlation structure and defining
     # correlation operator
-    p_spec = (lambda k: (42 / (k + 1) ** 3))
+    p_spec = (lambda k: (42/(k+1)**3))
 
     S = ift.create_power_operator(h_space, power_spectrum=p_spec)
 
@@ -35,7 +35,7 @@ if __name__ == "__main__":
     Instrument = ift.DiagonalOperator(diag)
 
     # Adding a harmonic transformation to the instrument
-    R = ift.ComposedOperator([fft, Instrument])
+    R = Instrument*fft
     signal_to_noise = 1.
     ndiag = ift.Field.full(s_space, ss.var()/signal_to_noise)
     N = ift.DiagonalOperator(ndiag.weight(1))
diff --git a/nifty/library/wiener_filter_curvature.py b/nifty/library/wiener_filter_curvature.py
index b6d72dfe1ff276f2e63b2f7141331c6be29204c3..cabe0ff9552f771b9f98ee4330cef08897362e1c 100644
--- a/nifty/library/wiener_filter_curvature.py
+++ b/nifty/library/wiener_filter_curvature.py
@@ -31,7 +31,7 @@ class WienerFilterCurvature(EndomorphicOperator):
 
     @property
     def domain(self):
-        return self.S.domain
+        return self._op.domain
 
     @property
     def capability(self):
diff --git a/nifty/operators/__init__.py b/nifty/operators/__init__.py
index b6d015c96410cb11847c0b6c3333181305543f9f..5c722d2bba208ea930e00a28a837b7e87bf4137a 100644
--- a/nifty/operators/__init__.py
+++ b/nifty/operators/__init__.py
@@ -6,7 +6,6 @@ from .fft_smoothing_operator import FFTSmoothingOperator
 from .direct_smoothing_operator import DirectSmoothingOperator
 from .fft_operator import FFTOperator
 from .inversion_enabler import InversionEnabler
-from .composed_operator import ComposedOperator
 from .response_operator import ResponseOperator
 from .laplace_operator import LaplaceOperator
 from .smoothness_operator import SmoothnessOperator
diff --git a/nifty/operators/composed_operator.py b/nifty/operators/composed_operator.py
deleted file mode 100644
index cc6872382bc0d1a31e534e63b1d691731e0f939b..0000000000000000000000000000000000000000
--- a/nifty/operators/composed_operator.py
+++ /dev/null
@@ -1,78 +0,0 @@
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program.  If not, see <http://www.gnu.org/licenses/>.
-#
-# Copyright(C) 2013-2017 Max-Planck-Society
-#
-# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
-# and financially supported by the Studienstiftung des deutschen Volkes.
-
-from .linear_operator import LinearOperator
-
-
-class ComposedOperator(LinearOperator):
-    """ NIFTY class for composed operators.
-
-    The  NIFTY composed operator class combines multiple linear operators.
-
-    Parameters
-    ----------
-    operators : tuple of NIFTy Operators
-        The tuple of LinearOperators.
-
-    Attributes
-    ----------
-    domain : DomainTuple
-        The NIFTy.space in which the operator is defined.
-    target : DomainTuple
-        The NIFTy.space in which the outcome of the operator lives
-    """
-
-    def __init__(self, operators):
-        super(ComposedOperator, self).__init__()
-
-        self._operator_store = ()
-        old_op = None
-        self._capability = operators[0].capability
-        for op in operators:
-            if not isinstance(op, LinearOperator):
-                raise TypeError("The elements of the operator list must be"
-                                "instances of the LinearOperator base class")
-            if old_op is not None and op.domain != old_op.target:
-                raise ValueError("incompatible domains")
-            self._operator_store += (op,)
-            self._capability &= op.capability
-            old_op = op
-        if self._capability == 0:
-            raise ValueError("composed operator does not support any mode")
-
-    @property
-    def domain(self):
-        return self._operator_store[0].domain
-
-    @property
-    def target(self):
-        return self._operator_store[-1].target
-
-    @property
-    def capability(self):
-        return self._capability
-
-    def apply(self, x, mode):
-        self._check_mode(mode)
-        if mode == self.TIMES or mode == self.ADJOINT_INVERSE_TIMES:
-            for op in self._operator_store:
-                x = op.apply(x, mode)
-        else:
-            for op in reversed(self._operator_store):
-                x = op.apply(x, mode)
-        return x
diff --git a/nifty/operators/fft_operator.py b/nifty/operators/fft_operator.py
index 4ce98e3694c637255cd5041f00d14a5ef414ce9c..79d41a7f8ea3919614e0dabfa396950556bbe725 100644
--- a/nifty/operators/fft_operator.py
+++ b/nifty/operators/fft_operator.py
@@ -54,20 +54,6 @@ class FFTOperator(LinearOperator):
         For GLSpace, HPSpace, and LMSpace, a sensible (but not unique)
         co-domain is chosen that should work satisfactorily in most situations,
         but for full control, the user should explicitly specify a codomain.
-
-    Attributes
-    ----------
-    domain: Tuple of Spaces
-        The domain of the data that is input by "times" and output by
-        "adjoint_times".
-    target: Tuple of Spaces
-        The domain of the data that is output by "times" and input by
-        "adjoint_times".
-
-    Raises
-    ------
-    ValueError:
-        if "domain" or "target" are not of the proper type.
     """
 
     def __init__(self, domain, target=None, space=None):
@@ -96,7 +82,8 @@ class FFTOperator(LinearOperator):
         else:
             self._trafo = SphericalTransformation(pdom, hdom, self._space)
 
-    def _times_helper(self, x):
+    def apply(self, x, mode):
+        self._check_input(x, mode)
         if np.issubdtype(x.dtype, np.complexfloating):
             res = (self._trafo.transform(x.real) +
                    1j * self._trafo.transform(x.imag))
@@ -104,10 +91,6 @@ class FFTOperator(LinearOperator):
             res = self._trafo.transform(x)
         return res
 
-    def apply(self, x, mode):
-        self._check_input(x, mode)
-        return self._times_helper(x)
-
     @property
     def domain(self):
         return self._domain
diff --git a/nifty/operators/response_operator.py b/nifty/operators/response_operator.py
index a348fd17300e0451afa72bd5470bda573c99cdc8..f3af37923a2ae9c23119054d90329f295e5fa060 100644
--- a/nifty/operators/response_operator.py
+++ b/nifty/operators/response_operator.py
@@ -2,8 +2,9 @@ from builtins import range
 from .. import Field, FieldArray, DomainTuple
 from .linear_operator import LinearOperator
 from .fft_smoothing_operator import FFTSmoothingOperator
-from .composed_operator import ComposedOperator
 from .diagonal_operator import DiagonalOperator
+from .scaling_operator import ScalingOperator
+import numpy as np
 
 
 class ResponseOperator(LinearOperator):
@@ -23,47 +24,39 @@ class ResponseOperator(LinearOperator):
         Defines the smoothing length of the operator for each space it lives on
     exposure : list(np.float)
         Defines the exposure of the operator for each space it lives on
-    default_spaces : tuple of ints *optional*
-        Defines on which space(s) of a given field the Operator acts by
-        default (default: None)
-
-    Attributes
-    ----------
-    domain : DomainTuple
-        The domain on which the Operator's input Field lives.
-    target : DomainTuple
-        The domain in which the outcome of the operator lives.
-
-    Raises
-    ------
-    ValueError:
-        raised if:
-            * len of sigma-list and exposure-list are not equal
+    spaces : tuple of ints *optional*
+        Defines on which subdomain the individual components act
+        (default: None)
     """
 
-    def __init__(self, domain, sigma=[1.], exposure=[1.], spaces=None):
+    def __init__(self, domain, sigma, exposure, spaces=None):
         super(ResponseOperator, self).__init__()
 
-        if len(sigma) != len(exposure):
-            raise ValueError("Length of smoothing kernel and length of"
-                             "exposure do not match")
-        nsigma = len(sigma)
-
         self._domain = DomainTuple.make(domain)
 
         if spaces is None:
-            spaces = range(len(self._domain))
-
-        kernel_smoothing = [
-            FFTSmoothingOperator(self._domain, sigma[x], space=spaces[x])
-            for x in range(nsigma)]
-        self._composed_kernel = ComposedOperator(kernel_smoothing)
-
-        kernel_exposure = [
-            DiagonalOperator(Field(self._domain[spaces[x]], exposure[x]),
-                             domain=self._domain, spaces=(spaces[x],))
-            for x in range(nsigma)]
-        self._composed_exposure = ComposedOperator(kernel_exposure)
+            spaces = tuple(range(len(self._domain)))
+        if len(sigma) != len(exposure) or len(sigma) != len(spaces):
+            raise ValueError("length mismatch between sigma, exposure, "
+                             "and spaces ")
+        ncomp = len(sigma)
+        if ncomp == 0:
+            raise ValueError("Empty response operator not allowed")
+
+        self._kernel = None
+        for i in range(ncomp):
+            op = FFTSmoothingOperator(self._domain, sigma[i], space=spaces[i])
+            self._kernel = op if self._kernel is None else op*self._kernel
+
+        self._exposure = None
+        for i in range(ncomp):
+            if np.isscalar(exposure[i]):
+                op = ScalingOperator(exposure[i], self._domain)
+            elif isinstance(exposure[i], Field):
+                op = DiagonalOperator(exposure[i], domain=self._domain,
+                                      spaces=(spaces[i],))
+            self._exposure = op if self._exposure is None\
+                else op*self._exposure
 
         target_list = [FieldArray(self._domain[i].shape) for i in spaces]
         self._target = DomainTuple.make(target_list)
@@ -81,19 +74,17 @@ class ResponseOperator(LinearOperator):
         return self.TIMES | self.ADJOINT_TIMES
 
     def _times(self, x):
-        res = self._composed_kernel.times(x)
-        res = self._composed_exposure.times(res)
+        res = self._exposure(self._kernel(x))
         # removing geometric information
         return Field(self._target, val=res.val)
 
     def _adjoint_times(self, x):
         # setting correct spaces
         res = Field(self.domain, val=x.val)
-        res = self._composed_exposure.adjoint_times(res)
+        res = self._exposure.adjoint_times(res)
         res = res.weight(power=-1)
-        return self._composed_kernel.adjoint_times(res)
+        return self._kernel.adjoint_times(res)
 
     def apply(self, x, mode):
         self._check_input(x, mode)
         return self._times(x) if mode == self.TIMES else self._adjoint_times(x)
-
diff --git a/nifty/operators/smoothness_operator.py b/nifty/operators/smoothness_operator.py
index 4308427310b3d83d7f5db5fada664cb6c93a6220..4f307a7842c1533b32106ecbe8193c9ddf9cab5b 100644
--- a/nifty/operators/smoothness_operator.py
+++ b/nifty/operators/smoothness_operator.py
@@ -38,6 +38,7 @@ class SmoothnessOperator(EndomorphicOperator):
     def domain(self):
         return self._laplace._domain
 
+    # MR FIXME: shouldn't this operator actually be self-adjoint?
     @property
     def capability(self):
         return self.TIMES
diff --git a/nifty/sugar.py b/nifty/sugar.py
index 98613ec77b99c79f1f728843552591d46606de89..f38fbd98b54ee0faa6ddf138e6ef0b43604716fd 100644
--- a/nifty/sugar.py
+++ b/nifty/sugar.py
@@ -17,7 +17,7 @@
 # and financially supported by the Studienstiftung des deutschen Volkes.
 
 import numpy as np
-from . import Space, PowerSpace, Field, ComposedOperator, DiagonalOperator,\
+from . import Space, PowerSpace, Field, DiagonalOperator,\
               PowerProjectionOperator, FFTOperator, sqrt, DomainTuple, dobj,\
               utilities
 
@@ -241,11 +241,10 @@ def create_power_operator(domain, power_spectrum, space=None, dtype=None):
 
 
 def create_composed_fft_operator(domain, codomain=None, all_to='other'):
-    fft_op_list = []
-
     if codomain is None:
         codomain = [None]*len(domain)
     tdom = list(domain)
+    res = None
     for i, space in enumerate(domain):
         if not isinstance(space, Space):
             continue
@@ -256,6 +255,9 @@ def create_composed_fft_operator(domain, codomain=None, all_to='other'):
                 tgt = domain[i].get_default_codomain()
             else:
                 tgt = codomain[i]
-            fft_op_list += [FFTOperator(domain=tdom, target=tgt, space=i)]
+            op = FFTOperator(domain=tdom, target=tgt, space=i)
+            res = op if res is None else op*res
             tdom[i] = tgt
-    return ComposedOperator(fft_op_list)
+    if res is None:
+        raise ValueError("empty operator")
+    return res
diff --git a/test/test_operators/test_composed_operator.py b/test/test_operators/test_composed_operator.py
index 11bb419f37e383a49b87b612390520f2fc91a3c8..c82aa1ce684375b1b5e82e854f805ffdf6814c80 100644
--- a/test/test_operators/test_composed_operator.py
+++ b/test/test_operators/test_composed_operator.py
@@ -17,7 +17,7 @@ class ComposedOperator_Tests(unittest.TestCase):
         op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
         op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
 
-        op = ift.ComposedOperator((op1, op2))
+        op = op2*op1
 
         rand1 = ift.Field.from_random('normal', domain=(space1, space2))
         rand2 = ift.Field.from_random('normal', domain=(space1, space2))
@@ -34,7 +34,7 @@ class ComposedOperator_Tests(unittest.TestCase):
         op1 = ift.DiagonalOperator(diag1, cspace, spaces=(0,))
         op2 = ift.DiagonalOperator(diag2, cspace, spaces=(1,))
 
-        op = ift.ComposedOperator((op1, op2))
+        op = op2*op1
 
         rand1 = ift.Field.from_random('normal', domain=(space1, space2))
         tt1 = op.inverse_times(op.times(rand1))