diff --git a/src/itensor.jl b/src/itensor.jl
index 293acdc3ce..7c9e374f63 100644
--- a/src/itensor.jl
+++ b/src/itensor.jl
@@ -4,6 +4,7 @@ export ITensor,
        combinedindex,
        delta,
        δ,
+       exp,
        replaceindex!,
        inds,
        isNull,
@@ -424,6 +425,28 @@ end
 
 dot(A::ITensor,B::ITensor) = scalar(dag(A)*B)
 
+import LinearAlgebra.exp
+
+"""
+    exp(A::ITensor, Lis::IndexSet; hermitian = false)
+Compute the exponential of the tensor `A` by treating it as a matrix ``A_{lr}`` with
+the left index `l` running over all indices in `Lis` and `r` running over all
+indices not in `Lis`. Must have `dim(Lis) == dim(inds(A))/dim(Lis)` for the exponentiation to
+be defined.
+When `hermitian=true` the exponential of `Hermitian(reshape(A, dim(Lis), :))` is
+computed internally.
+"""
+function exp(A::ITensor, Lis::IndexSet; hermitian = false)
+  (dim(Lis) == dim(inds(A))/dim(Lis)) || throw(DimensionMismatch("dimension of the left index set `Lis` must be
+                                                                       equal to `dim(inds(A))/dim(Lis)`"))
+  A, Lis, Ris = _permute_for_factorize(A,Lis)
+  expAs = storage_exp(store(A), Lis,Ris, hermitian = hermitian)
+  return ITensor(inds(A),expAs)
+end
+
+
+expHermitian(A::ITensor,Lis::IndexSet) = exp(A,Lis, hermitian=true)
+
 #######################################################################
 #
 # In-place operations
diff --git a/src/storage/dense.jl b/src/storage/dense.jl
index 691fa21f0f..2f90c590eb 100644
--- a/src/storage/dense.jl
+++ b/src/storage/dense.jl
@@ -357,3 +357,8 @@ function storage_polar(Astore::Dense{T},
   return (Qis,Qstore,Pis,Pstore)
 end
 
+function storage_exp(As::Dense{T}, Lis,Ris; hermitian=false) where {T}
+  expAdata = ( hermitian ? Array(exp(Hermitian(reshape(data(As),dim(Lis),dim(Ris))))) :
+                           exp(reshape(data(As),dim(Lis),dim(Ris))) )
+  return Dense{T}(vec(expAdata))
+end
diff --git a/test/test_itensor.jl b/test/test_itensor.jl
index 481d912245..430b2e4e9a 100644
--- a/test/test_itensor.jl
+++ b/test/test_itensor.jl
@@ -172,6 +172,41 @@ end
   @test dot(A, B) == 11.0
 end
 
+@testset "exponentiate" begin
+  s1 = Index(2,"s1")
+  s2 = Index(2,"s2")
+  i1 = Index(2,"i1")
+  i2 = Index(2,"i2")
+  Amat = rand(2,2,2,2)
+  A = ITensor(Amat, i1,i2,s1,s2)
+
+  Aexp = exp(A,IndexSet(i1,i2))
+  Amatexp = reshape( exp(reshape(Amat,4,4)), 2,2,2,2)
+  Aexp_from_mat = ITensor(Amatexp, i1,i2,s1,s2)
+  @test Aexp ≈ Aexp_from_mat
+
+  #test that exponentiation works when indices need to be permuted
+  Aexp = exp(A,IndexSet(s1,s2))
+  Amatexp = Array( exp(  reshape(Amat,4,4))' )
+  Aexp_from_mat = ITensor(reshape(Amatexp,2,2,2,2), s1,s2,i1,i2)
+  @test Aexp ≈ Aexp_from_mat
+
+  #test exponentiation when hermitian=true is used
+  Amat = reshape(Amat, 4,4)
+  Amat = reshape( Amat + Amat' + randn(4,4)*1e-10 , 2,2,2,2)
+  A = ITensor(Amat, i1,i2,s1,s2)
+  Aexp = exp(A,IndexSet(i1,i2), hermitian=true)
+  Amatexp = Array(reshape( exp(Hermitian(reshape(Amat,4,4))), 2,2,2,2))
+  Aexp_from_mat = ITensor(Amatexp, i1,i2,s1,s2)
+  @test Aexp ≈ Aexp_from_mat
+
+
+
+  @test_throws DimensionMismatch exp(A,IndexSet(s1))
+
+end
+
+
 @testset "add and axpy" begin
   i = Index(2,"i")
   a = [1.0; 2.0]