Repository: dgasmith/opt_einsum

Issue:
Optimal path with constant arguments
I have an application where there is a large tensor contraction in the inner loop, with some arrays that are fixed and others that are updated between iterations. Operations only involving the constant terms could done ahead of time (outside the loop), but this may or may not actually improve performance.

It could be nice if opt-einsum had a way to mark some arguments as constant, in which case the cost of contracting them them could neglected when computing the cost of paths only involving these arguments.

Discussion:
codecov-io: # [Codecov](https://codecov.io/gh/dgasmith/opt_einsum/pull/30?src=pr&el=h1) Report
> Merging [#30](https://codecov.io/gh/dgasmith/opt_einsum/pull/30?src=pr&el=desc) into [master](https://codecov.io/gh/dgasmith/opt_einsum/commit/960cf7c45993077be6d4ad04a1564b6042f0e092?src=pr&el=desc) will **increase** coverage by `0.1%`.
> The diff coverage is `100%`.


dgasmith: Cool! I will think a bit about how to observe constant expressions that are "future" in the path. Standard graph techniques could do this trivially, but seems like a lot of code when we can probably just use a few loops for it.

Also in the first example I think the constants `kwarg` should be `constants=[2, 3, 4]`?
jcmgray: > Cool! I will think a bit about how to observe constant expressions that are "future" in the path. Standard graph techniques could do this trivially, but seems like a lot of code when we can probably just use a few loops for it.

Yeah I haven't thought too much about these 'future' constant operations yet, but suspect they are quite rare.

One thing that would make addressing this easier would be to store the tensors not in a ``list`` but in a ``dict`` so that they are identifiable by more than just their position (which changes, meaning its hard to swap the order of operations). This kind of thing might also e.g. allow caching in the ``'optimal'`` algorithm and also reduce the cost of updating potential contractions in ``greedy``. Quite a big backend change though! And probably not necessary for this PR.

> Also in the first example I think the constants ``kwarg`` should be ``constants=[2, 3, 4]``?

The supplied ops are ``shape, array, array, array, shape`` so ``[1, 2, 3]`` is correct (if we are using 0-indexing). A point of confusion might be that the ``repr`` has all the constants moved to the end? :

```
'ij,jk,kl,lm,mn->ni' -> 'ij,mn,[jk,kl,lm]->ni'
```
jcmgray: Yes think I'll leave future constants for moment, as like you say, a more graph like structure will make them easier to implement.

Here's some themes and variations on potential representations:
```python
# 1 - nothing
"ContractExpression('ij,jk,kl,lm,mn->ni', constants=(1, 2, 4))"

# 2 - individually wrapped
"ContractExpression('ij,[jk,kl],lm,[mn]->ni', constants=(1, 2, 4))"

# 3 - Moved to end
"ContractExpression('ij,lm[jk,kl,mn]->ni')"

# 4 - explicitly mentioned in contraction string #1
"ContractExpression('ij,lm->ni [Constant: jk,kl,mn]')"

# 5 - explicitly mentioned in contraction string #2
"ContractExpression('ij,lm,[Const:jk,kl,mn]->ni')"

# 6 - explicitly mentioned in contraction string #3
"ContractExpression('ij,lm,Const[jk,kl,mn]->ni')"
```
I'd favour 2 or maybe 5. Its certainly nice to keep the them in order, apart from making it less clear how many arguments the expression should now be called with (realistically if the user has created the expression they probably know though!).

> Since we are moving around the order we may have to do something like Input: A_ij,B_[jk],C_[kl],D_[lm],E_mn->ni. Need to think about this more.

My 2c on introducing new letters is that might be a bit confusing, especially as the indices themselves could potentially be any upper or lower case characters.

---

I've also made a few tweaks to the backend stuff (and also will submit another PR later with support for torch and eager-mode tensorflow) but so far the performance gains from compiling these expressions with constants is pretty good, especially with a GPU!


> On switch to a dictionary it is certainly something to think about it. It would be nice to remove all of the indexing that we have floating around. I think you already switched the greedy algorithm to a dynamic approach to reduce the cost, not sure how it would directly effect optimal. It is certainly a point that we could do some sort of memoization on the contraction finder.

Yeah the logic for a dict in the optimal case is just that all tensors and contraction costs are uniquely defined by their indices (up to permutations), so a lot of different results could be reused regardless of however factorially many paths reach a certain point.
dgasmith: I guess my main concern is what if someone writes the following:

```
>>> contract_expression("ij,ij,ij,ij->, *ops, constants=[1, 3])
"ContractExpression('ij,ij,[Const:ij,ij]->')"
```

It is an edge case, but I think the loss of information there is pretty rough. All things considered I would vote for #2, I believe it is quite a bit more elegant than my try.

--

Great, I had been wondering about JIT for these expressions for a bit now. One thing to consider is `numba` or `numexpr` as well. 
jcmgray: Yes that is a good edge case to point out! I've gone with number 2 as it does seem to be a good tradeoff.

I've also had a first run at some documentation. Let me know if anything could be clearer.

Was also thinking about adding this and ``opt_einsum`` 's other features to a kind of 'highlights' list on the main readme. Thoughts?  

----
> Great, I had been wondering about JIT for these expressions for a bit now. One thing to consider is numba or numexpr as well.

Yes - a recent ``numba`` release actually added support for proper n-dimensional reshaping and transposing of arrays which should mean a numba tensordot implementation is possible.
dgasmith: Adding to the readme and docs landing page would be good. However, I do worry a bit about overall length and the amount of content that we throw at people. I think these pages could be ~20-30% larger without too much issue. More than that is fine, but we would need to cut down on basic `opt_einsum` info in that case to keep the overall length rather short.

Overall this is looking great and is pretty close to being ready to go IMO.

Patch:
diff --git a/docs/source/backends.rst b/docs/source/backends.rst
index 11efc96..9f05977 100644
--- a/docs/source/backends.rst
+++ b/docs/source/backends.rst
@@ -137,7 +137,7 @@ If ``theano`` is installed, using it as backend is as simple as specifiying
     >>> shapes = (3, 200), (200, 300), (300, 4)
     >>> expr = oe.contract_expression("ab,bc,cd", *shapes)
     >>> expr
-    ContractExpression('ab,bc,cd')
+    <ContractExpression('ab,bc,cd')>
 
     >>> import numpy as np
     >>> # GPU advantage mainly for low precision numbers
diff --git a/docs/source/reusing_paths.rst b/docs/source/reusing_paths.rst
index 27a37ff..d920d4a 100644
--- a/docs/source/reusing_paths.rst
+++ b/docs/source/reusing_paths.rst
@@ -8,7 +8,7 @@ If you expect to repeatedly use a particular contraction it can make things simp
 
     >>> my_expr = oe.contract_expression("abc,cd,dbe->ea", (2, 3, 4), (4, 5), (5, 3, 6))
     >>> print(my_expr)
-    <ContractExpression> for 'abc,cd,dbe->ea':
+    <ContractExpression('abc,cd,dbe->ea')>
       1.  'dbe,cd->bce' [GEMM]
       2.  'bce,abc->ea' [GEMM]
 
@@ -26,3 +26,109 @@ The ``ContractExpression`` can be called with 3 arrays that match the original s
            [ 3.67772272,  5.46727192]])
 
 Note that few checks are performed when calling the expression, and while it will work for a set of arrays with the same ranks as the original shapes but differing sizes, it might no longer be optimal.
+
+
+====================
+Specifying Constants
+====================
+
+Often one generates contraction expressions where some of the tensor arguments
+will remain *constant* across many calls.
+:func:`~opt_einsum.contract_expression` allows you to specify the indices of
+these constant arguments, allowing ``opt_einsum`` to build and then reuse as
+many constant contractions as possible. Take for example the equation:
+
+.. code:: python
+
+    >>> eq = "ij,jk,kl,lm,mn->ni"
+
+where we know that *only* the first and last tensors will vary between calls.
+We can specify this by marking the middle three as constant - we then need to
+supply the actual arrays rather than just the shapes to
+:func:`~opt_einsum.contract_expression`:
+
+.. code:: python
+
+    >>> #           A       B       C       D       E
+    >>> shapes = [(9, 5), (5, 5), (5, 5), (5, 5), (5, 8)]
+
+    >>> # mark the middle three arrays as constant
+    >>> constants = [1, 2, 3]
+
+    >>> # generate the constant arrays
+    >>> B, C, D = [np.random.randn(*shapes[i]) for i in constants]
+
+    >>> # supplied ops are now mix of shapes and arrays
+    >>> ops = (9, 5), B, C, D, (5, 8)
+
+    >>> expr = oe.contract_expression(eq, *ops, constants=constants)
+    >>> expr
+    <ContractExpression('ij,[jk,kl,lm],mn->ni', constants=[1, 2, 3])>
+
+The expression now only takes the remaining two arrays as arguments (the
+tensors with ``'ij'`` and ``'mn'`` indices), and will store as many resuable
+constant contractions as possible.
+
+.. code:: python
+
+    >>> A1, E1 = np.random.rand(*shapes[0]), np.random.rand(*shapes[-1])
+    >>> out1 = expr(A1, E1)
+    >>> out1.shap
+    (8, 9)
+
+    >>> A2, E2 = np.random.rand(*shapes[0]), np.random.rand(*shapes[-1])
+    >>> out2 = expr(A2, E2)
+    >>> out2.shape
+    (8, 9)
+
+    >>> np.allclose(out1, out2)
+    False
+
+    >>> print(expr)
+    <ContractExpression('ij,[jk,kl,lm],mn->ni', constants=[1, 2, 3])>
+      1.  'jm,mn->jn' [GEMM]
+      2.  'jn,ij->ni' [GEMM]
+
+Where we can see that the expression now only has to perform
+two contractions to compute the output.
+
+.. note::
+
+    The constant part of an expression is lazily generated upon first call,
+    (specific a particular backend) though it can be explicitly built with call
+    to :meth:`~opt_einsum.contract.ContractExpression.evaluate_constants`.
+
+We can confirm the advantage of using expressions and constants by timing the
+following scenarios, first setting
+``A = np.random.rand(*shapes[0])`` and ``E = np.random.rand(*shapes[-1])``.
+
+- **contract from scratch:**
+
+.. code:: python
+
+    >>> %timeit oe.contract(eq, A, B, C, D, E)
+    239 µs ± 5.06 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
+
+- **contraction with an expression but no constants:**
+
+.. code:: python
+
+    >>> expr_no_consts = oe.contract_expression(eq, *shapes)
+    >>> %timeit expr_no_consts(A, B, C, D, E)
+    76.7 µs ± 2.47 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
+
+- **contraction with an expression and constants marked:**
+
+.. code:: python
+
+    >>> %timeit expr(A, E)
+    40.8 µs ± 1.22 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
+
+Although this gives us a rough idea, of course the efficiency savings are
+hugely dependent on the size of the contraction and number of possible constant
+contractions.
+
+We also note that even if there are *no* constant contractions to perform, it
+can be very advantageous to specify constant tensors for particular backends.
+For instance, if a GPU backend is used, the constant tensors will be kept on
+the device rather than being transfered each time.
diff --git a/opt_einsum/backends.py b/opt_einsum/backends.py
index 2a08371..e133d6a 100644
--- a/opt_einsum/backends.py
+++ b/opt_einsum/backends.py
@@ -62,72 +62,130 @@ def has_einsum(backend):
 
 # Tensorflow
 
-def convert_arrays_to_tensorflow(arrays):
-    """Convert numpy arrays to ``tensorflow.placeholder`` instances.
+def to_tensorflow(array, constant=False):
+    """Convert a numpy array to a ``tensorflow.placeholder`` instance.
     """
     import tensorflow
-    return [tensorflow.placeholder(x.dtype, x.shape) for x in arrays]
+
+    if isinstance(array, numpy.ndarray):
+        if constant:
+            return tensorflow.constant(array, array.dtype, array.shape)
+
+        return tensorflow.placeholder(array.dtype, array.shape)
+
+    return array
 
 
 def build_tensorflow_expression(arrays, expr):
     """Build a tensorflow function based on ``arrays`` and ``expr``.
     """
     import tensorflow
-    placeholders = convert_arrays_to_tensorflow(arrays)
-    graph = expr._normal_contract(placeholders, backend='tensorflow')
+
+    placeholders = [to_tensorflow(array) for array in arrays]
+    graph = expr._contract(placeholders, backend='tensorflow')
 
     def tensorflow_contract(*arrays):
         session = tensorflow.get_default_session()
-        return session.run(graph, feed_dict=dict(zip(placeholders, arrays)))
+        # only want to feed placeholders - constant tensors already have values
+        feed_dict = {p: a for p, a in zip(placeholders, arrays) if p.op.type == 'Placeholder'}
+        return session.run(graph, feed_dict=feed_dict)
 
     return tensorflow_contract
 
 
+def evaluate_constants_tensorflow(const_arrays, expr):
+    """Convert constant arguments to tensorflow constants, and perform any
+    possible constant contractions. Requires evaluating a tensorflow graph.
+    """
+    import tensorflow
+
+    # compute the partial graph of new inputs
+    const_arrays = [to_tensorflow(x, constant=True) for x in const_arrays]
+    new_ops, new_contraction_list = expr(*const_arrays, backend='tensorflow', evaluate_constants=True)
+
+    # evaluate the new inputs and convert to tensorflow constants
+    session = tensorflow.get_default_session()
+    new_ops = [None if x is None else to_tensorflow(session.run(x), constant=True) for x in new_ops]
+
+    return new_ops, new_contraction_list
+
+
 # Theano
 
-def convert_arrays_to_theano(arrays):
-    """Convert numpy arrays to ``theano.tensor.TensorType`` instances.
+def to_theano(array, constant=False):
+    """Convert a numpy array to ``theano.tensor.TensorType`` instance.
     """
     import theano
-    return [theano.tensor.TensorType(dtype=x.dtype, broadcastable=[False] * len(x.shape))() for x in arrays]
+
+    if isinstance(array, numpy.ndarray):
+        if constant:
+            return theano.tensor.constant(array)
+
+        return theano.tensor.TensorType(dtype=array.dtype, broadcastable=[False] * len(array.shape))()
+
+    return array
 
 
 def build_theano_expression(arrays, expr):
     """Build a theano function based on ``arrays`` and ``expr``.
     """
     import theano
-    in_vars = convert_arrays_to_theano(arrays)
-    out_var = expr._normal_contract(in_vars, backend='theano')
-    graph = theano.function(in_vars, out_var)
+
+    in_vars = [to_theano(array) for array in arrays]
+    out_var = expr._contract(in_vars, backend='theano')
+
+    # don't supply constants to graph
+    graph_ins = [x for x in in_vars if not isinstance(x, theano.tensor.TensorConstant)]
+    graph = theano.function(graph_ins, out_var)
 
     def theano_contract(*arrays):
-        return graph(*arrays)
+        return graph(*[x for x in arrays if not isinstance(x, theano.tensor.TensorConstant)])
 
     return theano_contract
 
 
+def evaluate_constants_theano(const_arrays, expr):
+    # compute the partial graph of new inputs
+    const_arrays = [to_theano(x, constant=True) for x in const_arrays]
+    new_ops, new_contraction_list = expr(*const_arrays, backend='theano', evaluate_constants=True)
+
+    # evaluate the new inputs and convert to theano shared tensors
+    new_ops = [None if x is None else to_theano(x.eval(), constant=True) for x in new_ops]
+
+    return new_ops, new_contraction_list
+
+
 # Cupy
 
-def convert_arrays_to_cupy(arrays):  # pragma: no cover
-    """Convert numpy arrays to ``cupy.ndarray`` instances.
-    """
+def to_cupy(array):  # pragma: no cover
     import cupy
-    return [cupy.asarray(x) for x in arrays]
+
+    if isinstance(array, numpy.ndarray):
+        return cupy.asarray(array)
+
+    return array
 
 
 def build_cupy_expression(_, expr):  # pragma: no cover
     """Build a cupy function based on ``arrays`` and ``expr``.
     """
-    import cupy
 
     def cupy_contract(*arrays):
-        cupy_arrays = convert_arrays_to_cupy(arrays)
-        cupy_out = expr._normal_contract(cupy_arrays, backend='cupy')
-        return cupy.asnumpy(cupy_out)
+        cupy_arrays = [to_cupy(x) for x in arrays]
+        cupy_out = expr._contract(cupy_arrays, backend='cupy')
+        return cupy_out.get()
 
     return cupy_contract
 
 
+def evaluate_constants_cupy(const_arrays, expr):  # pragma: no cover
+    """Convert constant arguments to cupy arrays, and perform any possible
+    constant contractions.
+    """
+    const_arrays = [to_cupy(x) for x in const_arrays]
+    return expr(*const_arrays, backend='cupy', evaluate_constants=True)
+
+
 # Dispatch to correct expression backend
 #    these are the backends which support explicit to-and-from numpy conversion
 CONVERT_BACKENDS = {
@@ -137,8 +195,22 @@ CONVERT_BACKENDS = {
 }
 
 
+PARSE_CONSTS_BACKENDS = {
+    'tensorflow': evaluate_constants_tensorflow,
+    'theano': evaluate_constants_theano,
+    'cupy': evaluate_constants_cupy,
+}
+
+
 def build_expression(backend, arrays, expr):
     """Build an expression, based on ``expr`` and initial arrays ``arrays``,
     that evaluates using backend ``backend``.
     """
     return CONVERT_BACKENDS[backend](arrays, expr)
+
+
+def evaluate_constants(backend, arrays, expr):
+    """Convert constant arrays to the correct backend, and perform as much of
+    the contraction of ``expr`` with these as possible.
+    """
+    return PARSE_CONSTS_BACKENDS[backend](arrays, expr)
diff --git a/opt_einsum/contract.py b/opt_einsum/contract.py
index 651600e..e1fc026 100644
--- a/opt_einsum/contract.py
+++ b/opt_einsum/contract.py
@@ -402,7 +402,8 @@ def contract(*operands, **kwargs):
     use_blas = kwargs.pop('use_blas', True)
     memory_limit = kwargs.pop('memory_limit', None)
     backend = kwargs.pop('backend', 'numpy')
-    gen_expression = kwargs.pop('gen_expression', False)
+    gen_expression = kwargs.pop('_gen_expression', False)
+    constants_dict = kwargs.pop('_constants_dict', {})
 
     # Make sure remaining keywords are valid for einsum
     unknown_kwargs = [k for (k, v) in kwargs.items() if k not in valid_einsum_kwargs]
@@ -418,17 +419,17 @@ def contract(*operands, **kwargs):
 
     # check if performing contraction or just building expression
     if gen_expression:
-        return ContractExpression(full_str, contraction_list, **einsum_kwargs)
+        return ContractExpression(full_str, contraction_list, constants_dict, **einsum_kwargs)
 
     return _core_contract(operands, contraction_list, backend=backend, **einsum_kwargs)
 
 
-def _core_contract(operands, contraction_list, backend='numpy', **einsum_kwargs):
+def _core_contract(operands, contraction_list, backend='numpy', evaluate_constants=False, **einsum_kwargs):
     """Inner loop used to perform an actual contraction given the output
     from a ``contract_path(..., einsum_call=True)`` call.
     """
 
-    # Special handeling if out is specified
+    # Special handling if out is specified
     out_array = einsum_kwargs.pop('out', None)
     specified_out = out_array is not None
 
@@ -438,9 +439,13 @@ def _core_contract(operands, contraction_list, backend='numpy', **einsum_kwargs)
     # Start contraction loop
     for num, contraction in enumerate(contraction_list):
         inds, idx_rm, einsum_str, remaining, blas_flag = contraction
-        tmp_operands = []
-        for x in inds:
-            tmp_operands.append(operands.pop(x))
+
+        # check if we are performing the pre-pass of an expression with constants,
+        #     if so, break out upon finding first non-constant (None) operand
+        if evaluate_constants and any(operands[x] is None for x in inds):
+            return operands, contraction_list[num:]
+
+        tmp_operands = [operands.pop(x) for x in inds]
 
         # Do we need to deal with the output?
         handle_out = specified_out and ((num + 1) == len(contraction_list))
@@ -452,9 +457,7 @@ def _core_contract(operands, contraction_list, backend='numpy', **einsum_kwargs)
             input_str, results_index = einsum_str.split('->')
             input_left, input_right = input_str.split(',')
 
-            tensor_result = input_left + input_right
-            for s in idx_rm:
-                tensor_result = tensor_result.replace(s, "")
+            tensor_result = "".join(s for s in input_left + input_right if s not in idx_rm)
 
             # Find indices to contract over
             left_pos, right_pos = [], []
@@ -483,7 +486,7 @@ def _core_contract(operands, contraction_list, backend='numpy', **einsum_kwargs)
             # Do the contraction
             new_view = _einsum(einsum_str, *tmp_operands, backend=backend, **einsum_kwargs)
 
-        # Append new items and derefernce what we can
+        # Append new items and dereference what we can
         operands.append(new_view)
         del tmp_operands, new_view
 
@@ -493,35 +496,109 @@ def _core_contract(operands, contraction_list, backend='numpy', **einsum_kwargs)
         return operands[0]
 
 
+def format_const_einsum_str(einsum_str, constants):
+    """Add brackets to the constant terms in ``einsum_str``. For example:
+
+        >>> format_const_einsum_str('ab,bc,cd->ad', [0, 2])
+        'bc,[ab,cd]->ad'
+
+    No-op if there are no constants.
+    """
+    if not constants:
+        return einsum_str
+
+    if "->" in einsum_str:
+        lhs, rhs = einsum_str.split('->')
+        arrow = "->"
+    else:
+        lhs, rhs, arrow = einsum_str, "", ""
+
+    wrapped_terms = ["[{}]".format(t) if i in constants else t for i, t in enumerate(lhs.split(','))]
+
+    formatted_einsum_str = "{}{}{}".format(','.join(wrapped_terms), arrow, rhs)
+
+    # merge adjacent constants
+    formatted_einsum_str = formatted_einsum_str.replace("],[", ',')
+    return formatted_einsum_str
+
+
 class ContractExpression:
     """Helper class for storing an explicit ``contraction_list`` which can
     then be repeatedly called solely with the array arguments.
     """
 
-    def __init__(self, contraction, contraction_list, **einsum_kwargs):
-        self.contraction = contraction
+    def __init__(self, contraction, contraction_list, constants_dict, **einsum_kwargs):
         self.contraction_list = contraction_list
         self.einsum_kwargs = einsum_kwargs
-        self.num_args = len(contraction.split('->')[0].split(','))
+        self.contraction = format_const_einsum_str(contraction, constants_dict.keys())
+
+        # need to know _full_num_args to parse constants with, and num_args to call with
+        self._full_num_args = contraction.count(',') + 1
+        self.num_args = self._full_num_args - len(constants_dict)
+
+        # likewise need to know full contraction list
+        self._full_contraction_list = contraction_list
 
-    def _normal_contract(self, arrays, out=None, backend='numpy'):
+        self._constants_dict = constants_dict
+        self._evaluated_constants = {}
+        self._backend_expressions = {}
+
+    def evaluate_constants(self, backend='numpy'):
+        """Convert any constant operands to the correct backend form, and
+        perform as many contractions as possible to create a new list of
+        operands, stored in ``self._evaluated_constants[backend]``. This also
+        makes sure ``self.contraction_list`` only contains the remaining,
+        non-const operations.
+        """
+        # prepare a list of operands, with `None` for non-consts
+        tmp_const_ops = [self._constants_dict.get(i, None) for i in range(self._full_num_args)]
+
+        # get the new list of operands with constant operations performed, and remaining contractions
+        new_ops, new_contraction_list = self(*tmp_const_ops, backend=backend, evaluate_constants=True)
+        self._evaluated_constants[backend] = new_ops
+        self.contraction_list = new_contraction_list
+
+    def _get_evaluated_constants(self, backend):
+        """Retrieve or generate the cached list of constant operators (mixed
+        in with None representing non-consts) and the remaining contraction
+        list.
+        """
+        try:
+            return self._evaluated_constants[backend]
+        except KeyError:
+            self.evaluate_constants(backend)
+            return self._evaluated_constants[backend]
+
+    def _get_backend_expression(self, arrays, backend):
+        try:
+            return self._backend_expressions[backend]
+        except KeyError:
+            fn = backends.build_expression(backend, arrays, self)
+            self._backend_expressions[backend] = fn
+            return fn
+
+    def _contract(self, arrays, out=None, backend='numpy', evaluate_constants=False):
         """The normal, core contraction.
         """
-        return _core_contract(list(arrays), self.contraction_list, out=out,
-                              backend=backend, **self.einsum_kwargs)
+        contraction_list = self._full_contraction_list if evaluate_constants else self.contraction_list
+
+        return _core_contract(list(arrays), contraction_list, out=out, backend=backend,
+                              evaluate_constants=evaluate_constants, **self.einsum_kwargs)
 
-    def _convert_contract(self, arrays, out, backend):
+    def _contract_with_conversion(self, arrays, out, backend, evaluate_constants=False):
         """Special contraction, i.e. contraction with a different backend
-        but converting to and from that backend. Checks for
-        ``self._{backend}_contract``, generates it if is missing, then calls it
+        but converting to and from that backend. Retrieves or generates a
+        cached expression using ``arrays`` as templates, then calls it
         with ``arrays``.
-        """
-        convert_fn = "_{}_contract".format(backend)
 
-        if not hasattr(self, convert_fn):
-            setattr(self, convert_fn, backends.build_expression(backend, arrays, self))
+        If ``evaluate_constants=True``, perform a partial contraction that
+        prepares the constant tensors and operations with the right backend.
+        """
+        # convert consts to correct type & find reduced contraction list
+        if evaluate_constants:
+            return backends.evaluate_constants(backend, arrays, self)
 
-        result = getattr(self, convert_fn)(*arrays)
+        result = self._get_backend_expression(arrays, backend)(*arrays)
 
         if out is not None:
             out[()] = result
@@ -543,24 +620,34 @@ class ContractExpression:
             are supplied then try to convert them to and from the correct
             backend array type.
         """
-
-        if len(arrays) != self.num_args:
-            raise ValueError("This `ContractExpression` takes exactly %s array arguments "
-                             "but received %s." % (self.num_args, len(arrays)))
-
-        backend = kwargs.pop('backend', 'numpy')
         out = kwargs.pop('out', None)
+        backend = kwargs.pop('backend', 'numpy')
+        evaluate_constants = kwargs.pop('evaluate_constants', False)
+
         if kwargs:
             raise ValueError("The only valid keyword arguments to a `ContractExpression` "
                              "call are `out=` or `backend=`. Got: %s." % kwargs)
 
+        correct_num_args = self._full_num_args if evaluate_constants else self.num_args
+
+        if len(arrays) != correct_num_args:
+            raise ValueError("This `ContractExpression` takes exactly %s array arguments "
+                             "but received %s." % (self.num_args, len(arrays)))
+
+        if self._constants_dict and not evaluate_constants:
+            # fill in the missing non-constant terms with newly supplied arrays
+            ops_var, ops_const = iter(arrays), self._get_evaluated_constants(backend)
+            ops = [next(ops_var) if op is None else op for op in ops_const]
+        else:
+            ops = arrays
+
         try:
             # Check if the backend requires special preparation / calling
             #   but also ignore non-numpy arrays -> assume user wants same type back
-            if backend in backends.CONVERT_BACKENDS and isinstance(arrays[0], np.ndarray):
-                return self._convert_contract(arrays, out, backend)
+            if backend in backends.CONVERT_BACKENDS and any(isinstance(x, np.ndarray) for x in arrays):
+                return self._contract_with_conversion(ops, out, backend, evaluate_constants=evaluate_constants)
 
-            return self._normal_contract(arrays, out, backend)
+            return self._contract(ops, out, backend, evaluate_constants=evaluate_constants)
 
         except ValueError as err:
             original_msg = str(err.args) if err.args else ""
@@ -571,10 +658,14 @@ class ContractExpression:
             raise
 
     def __repr__(self):
-        return "ContractExpression('%s')" % self.contraction
+        if self._constants_dict:
+            constants_repr = ", constants={}".format(sorted(self._constants_dict))
+        else:
+            constants_repr = ""
+        return "<ContractExpression('{}'{})>".format(self.contraction, constants_repr)
 
     def __str__(self):
-        s = "<ContractExpression> for '%s':" % self.contraction
+        s = self.__repr__()
         for i, c in enumerate(self.contraction_list):
             s += "\n  %i.  " % (i + 1)
             s += "'%s'" % c[2] + (" [%s]" % c[-1] if c[-1] else "")
@@ -602,6 +693,13 @@ def contract_expression(subscripts, *shapes, **kwargs):
         Specifies the subscripts for summation.
     shapes : sequence of integer tuples
         Shapes of the arrays to optimize the contraction for.
+    constants : sequence of int, optional
+        The indices of any constant arguments in ``shapes``, in which case the
+        actual array should be supplied at that position rather than just a
+        shape. If these are specified, then constant parts of the contraction
+        between calls will be reused. Additionally, if a gpu-enabled backend is
+        used for example, then the constant tensors will be kept on the gpu,
+        minimizing transfers.
     kwargs :
         Passed on to ``contract_path`` or ``einsum``. See ``contract``.
 
@@ -621,15 +719,29 @@ def contract_expression(subscripts, *shapes, **kwargs):
     - The generated expression will work with any arrays which have
       the same rank (number of dimensions) as the original shapes, however, if
       the actual sizes are different, the expression may no longer be optimal.
+    - Constant operations will be computed upon first call with a particular
+      backend, then subsequently reused.
 
     Examples
     --------
 
-    >>> expr = contract_expression("ab,bc->ac", (3, 4), (4, 5))
-    >>> a, b = np.random.rand(3, 4), np.random.rand(4, 5)
-    >>> c = expr(a, b)
-    >>> np.allclose(c, a @ b)
-    True
+    Basic usage:
+
+        >>> expr = contract_expression("ab,bc->ac", (3, 4), (4, 5))
+        >>> a, b = np.random.rand(3, 4), np.random.rand(4, 5)
+        >>> c = expr(a, b)
+        >>> np.allclose(c, a @ b)
+        True
+
+    Supply ``a`` as a constant:
+
+        >>> expr = contract_expression("ab,bc->ac", a, (4, 5), constants=[0])
+        >>> expr
+        <ContractExpression('[ab],bc->ac', constants=[0])>
+
+        >>> c = expr(b)
+        >>> np.allclose(c, a @ b)
+        True
 
     """
     if not kwargs.get('optimize', True):
@@ -640,6 +752,14 @@ def contract_expression(subscripts, *shapes, **kwargs):
             raise ValueError("'%s' should only be specified when calling a "
                              "`ContractExpression`, not when building it." % arg)
 
-    dummy_arrays = [_ShapeOnly(s) for s in shapes]
+    kwargs['_gen_expression'] = True
+
+    # build dict of constant indices mapped to arrays
+    constants = kwargs.pop('constants', ())
+    constants_dict = {i: shapes[i] for i in constants}
+    kwargs['_constants_dict'] = constants_dict
+
+    # apart from constant arguments, make dummy arrays
+    dummy_arrays = [s if i in constants else _ShapeOnly(s) for i, s in enumerate(shapes)]
 
-    return contract(subscripts, *dummy_arrays, gen_expression=True, **kwargs)
+    return contract(subscripts, *dummy_arrays, **kwargs)


Test patch:
diff --git a/opt_einsum/tests/test_backends.py b/opt_einsum/tests/test_backends.py
index b2e43ba..19d72c2 100644
--- a/opt_einsum/tests/test_backends.py
+++ b/opt_einsum/tests/test_backends.py
@@ -6,6 +6,7 @@ from opt_einsum import contract, helpers, contract_expression, backends
 try:
     import tensorflow as tf
     found_tensorflow = True
+    sess = tf.Session()
 except ImportError:
     found_tensorflow = False
 
@@ -57,52 +58,131 @@ def test_tensorflow(string):
     shps = [v.shape for v in views]
     expr = contract_expression(string, *shps, optimize=True)
 
-    sess = tf.Session()
     with sess.as_default():
         expr(*views, backend='tensorflow', out=opt)
 
     assert np.allclose(ein, opt)
 
     # test non-conversion mode
-    tensorflow_views = backends.convert_arrays_to_tensorflow(views)
+    tensorflow_views = [backends.to_tensorflow(view) for view in views]
     expr(*tensorflow_views, backend='tensorflow')
 
 
-@pytest.mark.skipif(not found_cupy, reason="Cupy not installed.")
+@pytest.mark.skipif(not found_tensorflow, reason="Tensorflow not installed.")
+def test_tensorflow_with_constants():
+    eq = 'ij,jk,kl->li'
+    shapes = (2, 3), (3, 4), (4, 5)
+    constants = {0, 2}
+    ops = [np.random.rand(*shp) if i in constants else shp for i, shp in enumerate(shapes)]
+    var = np.random.rand(*shapes[1])
+
+    res_exp = contract(eq, ops[0], var, ops[2])
+
+    expr = contract_expression(eq, *ops, constants=constants)
+
+    # check tensorflow
+    with sess.as_default():
+        res_got = expr(var, backend='tensorflow')
+    assert 'tensorflow' in expr._evaluated_constants
+    assert np.allclose(res_exp, res_got)
+
+    # check can call with numpy still
+    res_got2 = expr(var, backend='numpy')
+    assert np.allclose(res_exp, res_got2)
+
+    # check tensorflow call returns tensorflow still
+    res_got3 = expr(backends.to_tensorflow(var), backend='tensorflow')
+    assert isinstance(res_got3, tf.Tensor)
+
+
+@pytest.mark.skipif(not found_theano, reason="Theano not installed.")
 @pytest.mark.parametrize("string", tests)
-def test_cupy(string):  # pragma: no cover
+def test_theano(string):
     views = helpers.build_views(string)
     ein = contract(string, *views, optimize=False, use_blas=False)
     shps = [v.shape for v in views]
 
     expr = contract_expression(string, *shps, optimize=True)
 
-    opt = expr(*views, backend='cupy')
+    opt = expr(*views, backend='theano')
     assert np.allclose(ein, opt)
 
     # test non-conversion mode
-    cupy_views = backends.convert_arrays_to_cupy(views)
-    cupy_opt = expr(*cupy_views, backend='cupy')
-    assert isinstance(cupy_opt, cupy.ndarray)
-    assert np.allclose(ein, cupy.asnumpy(cupy_opt))
+    theano_views = [backends.to_theano(view) for view in views]
+    theano_opt = expr(*theano_views, backend='theano')
+    assert isinstance(theano_opt, theano.tensor.TensorVariable)
 
 
-@pytest.mark.skipif(not found_theano, reason="Theano not installed.")
+@pytest.mark.skipif(not found_theano, reason="theano not installed.")
+def test_theano_with_constants():
+    eq = 'ij,jk,kl->li'
+    shapes = (2, 3), (3, 4), (4, 5)
+    constants = {0, 2}
+    ops = [np.random.rand(*shp) if i in constants else shp for i, shp in enumerate(shapes)]
+    var = np.random.rand(*shapes[1])
+
+    res_exp = contract(eq, ops[0], var, ops[2])
+
+    expr = contract_expression(eq, *ops, constants=constants)
+
+    # check theano
+    res_got = expr(var, backend='theano')
+    assert 'theano' in expr._evaluated_constants
+    assert np.allclose(res_exp, res_got)
+
+    # check can call with numpy still
+    res_got2 = expr(var, backend='numpy')
+    assert np.allclose(res_exp, res_got2)
+
+    # check theano call returns theano still
+    res_got3 = expr(backends.to_theano(var), backend='theano')
+    assert isinstance(res_got3, theano.tensor.TensorVariable)
+
+
+@pytest.mark.skipif(not found_cupy, reason="Cupy not installed.")
 @pytest.mark.parametrize("string", tests)
-def test_theano(string):
+def test_cupy(string):  # pragma: no cover
     views = helpers.build_views(string)
     ein = contract(string, *views, optimize=False, use_blas=False)
     shps = [v.shape for v in views]
 
     expr = contract_expression(string, *shps, optimize=True)
 
-    opt = expr(*views, backend='theano')
+    opt = expr(*views, backend='cupy')
     assert np.allclose(ein, opt)
 
     # test non-conversion mode
-    theano_views = backends.convert_arrays_to_theano(views)
-    theano_opt = expr(*theano_views, backend='theano')
-    assert isinstance(theano_opt, theano.tensor.TensorVariable)
+    cupy_views = [backends.to_cupy(view) for view in views]
+    cupy_opt = expr(*cupy_views, backend='cupy')
+    assert isinstance(cupy_opt, cupy.ndarray)
+    assert np.allclose(ein, cupy.asnumpy(cupy_opt))
+
+
+@pytest.mark.skipif(not found_cupy, reason="Cupy not installed.")
+def test_cupy_with_constants():
+    eq = 'ij,jk,kl->li'
+    shapes = (2, 3), (3, 4), (4, 5)
+    constants = {0, 2}
+    ops = [np.random.rand(*shp) if i in constants else shp for i, shp in enumerate(shapes)]
+    var = np.random.rand(*shapes[1])
+
+    res_exp = contract(eq, ops[0], var, ops[2])
+
+    expr = contract_expression(eq, *ops, constants=constants)
+
+    # check cupy
+    res_got = expr(var, backend='cupy')
+    assert 'cupy' in expr._evaluated_constants
+    assert np.allclose(res_exp, res_got)
+
+    # check can call with numpy still
+    res_got2 = expr(var, backend='numpy')
+    assert np.allclose(res_exp, res_got2)
+
+    # check cupy call returns cupy still
+    res_got3 = expr(cupy.asarray(var), backend='cupy')
+    assert isinstance(res_got3, cupy.ndarray)
+    assert np.allclose(res_exp, res_got3.get())
 
 
 @pytest.mark.skipif(not found_dask, reason="Dask not installed.")
diff --git a/opt_einsum/tests/test_contract.py b/opt_einsum/tests/test_contract.py
index e9c84ce..c439326 100644
--- a/opt_einsum/tests/test_contract.py
+++ b/opt_einsum/tests/test_contract.py
@@ -205,3 +205,31 @@ def test_contract_expressions(string, optimize, use_blas, out_spec):
     assert string in expr.__repr__()
     assert string in expr.__str__()
 
+
+@pytest.mark.parametrize("string,constants", [
+    ('hbc,bdef,cdkj,ji,ikeh,lfo', [1, 2, 3, 4]),
+    ('bdef,cdkj,ji,ikeh,hbc,lfo', [0, 1, 2, 3]),
+    ('hbc,bdef,cdkj,ji,ikeh,lfo', [1, 2, 3, 4]),
+    ('hbc,bdef,cdkj,ji,ikeh,lfo', [1, 2, 3, 4]),
+    ('ijab,acd,bce,df,ef->ji', [1, 2, 3, 4]),
+    ('ab,cd,ad,cb', [1, 3]),
+    ('ab,bc,cd', [0, 1]),
+])
+def test_contract_expression_with_constants(string, constants):
+    views = helpers.build_views(string)
+    expected = contract(string, *views, optimize=False, use_blas=False)
+
+    shapes = [view.shape for view in views]
+
+    expr_args = []
+    ctrc_args = []
+    for i, (shape, view) in enumerate(zip(shapes, views)):
+        if i in constants:
+            expr_args.append(view)
+        else:
+            expr_args.append(shape)
+            ctrc_args.append(view)
+
+    expr = contract_expression(string, *expr_args, constants=constants)
+    out = expr(*ctrc_args)
+    assert np.allclose(expected, out)