'Solve for z in (alpha + beta.wd.wd^T)z = wd . x + alpha z_est - gamma using gradient descent.
Output goes in z. temp is a matrix to store wd^Twd.'

def IterateForZ(self, train=False):

    input_layer = self.input_layer
    code_layer = self.code_layer
    epsilon = 0.01
    steps = 20
    z = code_layer.state
    wd = self.decoder.params['weight']
    hyp = code_layer.hyperparams
    alpha = hyp.sc_alpha
    beta = hyp.sc_beta
    gamma = hyp.sc_gamma
    temp = code_layer.m_by_m
    temp2 = code_layer.deriv
    temp3 = code_layer.temp3
    grad = code_layer.grad
    z_est = code_layer.approximator
    avg_models = (hyp.dropout and ((not hyp.dropout) or (not train)))
    cm.dot(wd, wd.T, target=temp)
    temp.mult(beta)
    if avg_models:
        temp.mult(((1.0 - hyp.dropout_prob) ** 2))
        temp.mult_diagonal((1.0 / (1.0 - hyp.dropout_prob)))
    temp.add_diagonal(alpha)
    z_est.mult(alpha, target=temp2)
    if avg_models:
        temp2.add_dot(wd, input_layer.state, mult=(beta * (1.0 - hyp.dropout_prob)))
    elif hyp.dropout:
        temp2.add_dot(wd, input_layer.state, mult=beta)
        temp2.mult(code_layer.mask)
    else:
        temp2.add_dot(wd, input_layer.state, mult=beta)
    z.assign(z_est)
    for i in range(steps):
        cm.dot(temp, z, target=grad)
        grad.subtract(temp2)
        z.sign(target=temp3)
        grad.add_mult(temp3, alpha=gamma)
        if (hyp.dropout and train):
            grad.mult(code_layer.mask)
        z.add_mult(grad, alpha=(- epsilon))
