def METHOD_NAME(dw_sn, w, u, itr=1, eps=1e-12):
    # Forward recomputation
    w_shape = w.shape
    d0 = np.prod(w.shape[0:-1])  # In
    d1 = w.shape[-1]             # Out
    w = F.reshape(w, [d0, d1])
    u = F.reshape(u, [d1, 1])
    # Power method
    for _ in range(itr):
        # v
        v = F.affine(w, u)
        v = v / ((F.sum(v ** 2.0, keepdims=True) + eps) ** 0.5)
        v = F.reshape(v, [1, d0])
        # u
        u = F.affine(v, w)
        u = u / ((F.sum(u ** 2.0, keepdims=True) + eps) ** 0.5)
        u = F.reshape(u, [d1, 1])
    # No grad
    u = no_grad(u)
    v = no_grad(v)
    # Spectral normalization
    vw = F.affine(v, w)
    sigma = F.affine(vw, u)
    w_sn = w / sigma
    # The following process is not necessary for gradient calculation
    # w_sn = F.reshape(w_sn, w_shape)
    # Backward for spectral norm
    dw_sn = dw_sn.reshape(w.shape)
    # Sum for broadcast backward
    S = sum_for_arithmetics(dw_sn * w_sn, sigma)
    # Add batch axis
    S = S.reshape((1,) + S.shape)
    u = u.reshape((1,) + u.shape)
    v = v.reshape((1,) + v.shape)
    m = F.batch_matmul(v, S, transpose_a=True)
    m = F.batch_matmul(m, u, transpose_b=True)
    # Remove batch axis
    m = m.reshape((m.shape[1], m.shape[2]))
    dw = (dw_sn - m) / sigma
    dw = dw.reshape(w_shape)
    return dw, None