def projector_splitting_add(Y, delta, debug=False):
    """Compute Y + delta via the projector splitting scheme.

    This function implements the projector splitting scheme (section 4.2 of [1]).

    The result is a TT-tensor with the TT-ranks equal to the TT-ranks of Y."""
    # Get rid of redundant ranks (they cause technical difficulties).
    delta = delta.round(eps=0)
    numDims = delta.d
    assert(numDims == Y.d)
    modeSize = delta.n
    assert(modeSize == Y.n).all()
    coresDelta = tt.tensor.to_list(delta)
    coresY = tt.tensor.to_list(Y)
    # rhs[dim] is an (delta.rank_dim * Y.rank_dim) x 1 vector
    rhs = [None] * (numDims+1)
    rhs[numDims] = np.ones([1, 1])
    # Right to left sweep to orthogonalize the cores and prepare the rhs.
    for dim in xrange(numDims-1, 0, -1):
        # Right to left orthogonalization of the Y cores.
        coresY = cores_orthogonalization_step(coresY, dim, left_to_right=False)
        r1, n, r2 = coresY[dim].shape

        # rhs computation.
        coreProd = np.tensordot(coresDelta[dim], coresY[dim], axes=(1, 1))
        coreProd = np.transpose(coreProd, (0, 2, 1, 3))
        coreProd = reshape(coreProd, (delta.r[dim]*r1, delta.r[dim+1]*r2))
        rhs[dim] = np.dot(coreProd, rhs[dim+1])
    if debug:
        assert(np.allclose(Y.full(), tt.tensor.from_list(coresY).full()))

    # lsh is an Y.rank_dim x delta.rank_dim matrix.
    lhs = np.ones([1, 1])
    # s is an Y.rank_dim x Y.rank_dim matrix.
    s = np.ones([1, 1])
    # Left to right projector splitting sweep.
    for dim in xrange(numDims):
        # Y^+ (formula 4.10)
        cc = coresDelta[dim].copy()
        r1, n, r2 = coresY[dim].shape
        cc = np.tensordot(lhs, cc, 1)
        rhs[dim+1] = reshape(rhs[dim+1], (delta.r[dim+1], r2))
        cc = reshape(cc, (-1, delta.r[dim+1]))
        cc = np.dot(cc, rhs[dim+1])
        if debug:
            first = np.kron(np.eye(modeSize[dim]), left(tt.tensor.from_list(coresY), dim-1).T)
            second = np.dot(first, unfolding(delta, dim))
            explicit = np.dot(second, right(tt.tensor.from_list(coresY), dim+1))
            assert(np.allclose(explicit, cc))
        cc += reshape(np.tensordot(s, coresY[dim], 1), (-1, Y.r[dim+1]))
        if dim < numDims-1:
            cc, rr = np.linalg.qr(cc)
            # TODO: do we need to use r1 = cc.shape[1] here????
        cc = reshape(cc, coresY[dim].shape)
        coresY[dim] = cc.copy()

        if dim < numDims-1:
            coreProd = np.tensordot(coresY[dim], coresDelta[dim], axes=(1, 1))
            coreProd = np.transpose(coreProd, (0, 2, 1, 3))
            coreProd = reshape(coreProd, (r1*delta.r[dim], r2*delta.r[dim+1]))
            lhs = reshape(lhs, (1, -1))
            lhs = np.dot(lhs, coreProd)
            lhs = reshape(lhs, (r2, delta.r[dim+1]))

        if dim < numDims-1:
            # Y^- (formula 4.7)
            s = rr - np.dot(lhs, rhs[dim+1])
            if debug:
                first = left(tt.tensor.from_list(coresY), dim).T
                second = np.dot(first, unfolding(delta, dim))
                explicit = np.dot(second, right(tt.tensor.from_list(coresY), dim+1))
                assert(np.allclose(explicit, np.dot(lhs, rhs[dim+1])))

    return tt.tensor.from_list(coresY)
Compute Y + delta via the projector splitting scheme.

    This function implements the projector splitting scheme (section 4.2 of [1]).

    The result is a TT-tensor with the TT-ranks equal to the TT-ranks of Y.