diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 12ff74c5..2a9c3830 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -147,14 +147,6 @@ class kern(parameterised):
         """
         K1 = self.copy()
         K2 = other.copy()
-        K1.unconstrain('')
-        K2.unconstrain('')
-
-        prev_ties = K1.tied_indices + [arr + K1.Nparam for arr in K2.tied_indices]
-        K1.untie_everything()
-        K2.untie_everything()
-
-        D = K1.D + K2.D
 
         newkernparts = [product_orthogonal(k1,k2) for k1, k2 in itertools.product(K1.parts,K2.parts)]
 
@@ -164,9 +156,14 @@ class kern(parameterised):
             s1[sl1], s2[sl2] = [True], [True]
             slices += [s1+s2]
 
-        newkern = kern(D, newkernparts, slices)
+        newkern = kern(K1.D + K2.D, newkernparts, slices)
+        newkern._follow_constrains(K1,K2)
+        
+        return newkern
 
-        # create the ties
+    def _follow_constrains(self,K1,K2):
+
+        # Build the array that allows to go from the initial indices of the param to the new ones
         K1_param = []
         n = 0
         for k1 in K1.parts:
@@ -180,19 +177,40 @@ class kern(parameterised):
         index_param = []
         for p1 in K1_param:
             for p2 in K2_param:
-                index_param += [0] + p1[1:] + p2[1:]
+                index_param += p1 + p2
         index_param = np.array(index_param)
 
+        # Get the ties and constrains of the kernels before the multiplication
+        prev_ties = K1.tied_indices + [arr + K1.Nparam for arr in K2.tied_indices]
+
+        prev_constr_pos = np.append(K1.constrained_positive_indices, K1.Nparam + K2.constrained_positive_indices)
+        prev_constr_neg = np.append(K1.constrained_negative_indices, K1.Nparam + K2.constrained_negative_indices)
+
+        prev_constr_fix = K1.constrained_fixed_indices + [arr + K1.Nparam for arr in K2.constrained_fixed_indices]
+        prev_constr_fix_values = K1.constrained_fixed_values + K2.constrained_fixed_values
+
+        prev_constr_bou = K1.constrained_bounded_indices + [arr + K1.Nparam for arr in K2.constrained_bounded_indices]
+        prev_constr_bou_low = K1.constrained_bounded_lowers + K2.constrained_bounded_lowers
+        prev_constr_bou_upp = K1.constrained_bounded_uppers + K2.constrained_bounded_uppers
+
         # follow the previous ties
         for arr in prev_ties:
             for j in arr:
                 index_param[np.where(index_param==j)[0]] = arr[0]
 
-        # tie
-        for i in np.unique(index_param)[1:]:
-            newkern.tie_param(np.where(index_param==i)[0])
-
-        return newkern
+        # ties and constrains
+        for i in range(K1.Nparam + K2.Nparam):
+            index = np.where(index_param==i)[0]
+            if index.size > 1:
+                self.tie_param(index)
+        for i in prev_constr_pos:
+            self.constrain_positive(np.where(index_param==i)[0])
+        for i in prev_constr_neg:
+            self.constrain_neg(np.where(index_param==i)[0])
+        for j, i in enumerate(prev_constr_fix):
+            self.constrain_fixed(np.where(index_param==i)[0],prev_constr_fix_values[j])
+        for j, i in enumerate(prev_constr_bou):
+            self.constrain_bounded(np.where(index_param==i)[0],prev_constr_bou_low[j],prev_constr_bou_upp[j])
 
     def _get_params(self):
         return np.hstack([p._get_params() for p in self.parts])
diff --git a/GPy/kern/product_orthogonal.py b/GPy/kern/product_orthogonal.py
index 26456d8c..46055885 100644
--- a/GPy/kern/product_orthogonal.py
+++ b/GPy/kern/product_orthogonal.py
@@ -4,7 +4,7 @@
 from kernpart import kernpart
 import numpy as np
 import hashlib
-from scipy import integrate
+#from scipy import integrate # This may not be necessary (Nicolas, 20th Feb)
 
 class product_orthogonal(kernpart):
     """
@@ -16,13 +16,12 @@ class product_orthogonal(kernpart):
 
     """
     def __init__(self,k1,k2):
-        assert k1._get_param_names()[0] == 'variance' and k2._get_param_names()[0] == 'variance', "Error: The multipication of kernels is only defined when the first parameters of the kernels to multiply is the variance."
         self.D = k1.D + k2.D
-        self.Nparam = k1.Nparam + k2.Nparam - 1
+        self.Nparam = k1.Nparam + k2.Nparam
         self.name = k1.name + '<times>' + k2.name
         self.k1 = k1
         self.k2 = k2
-        self._set_params(np.hstack((k1._get_params()[0]*k2._get_params()[0], k1._get_params()[1:],k2._get_params()[1:])))
+        self._set_params(np.hstack((k1._get_params(),k2._get_params())))
 
     def _get_params(self):
         """return the value of the parameters."""
@@ -30,14 +29,14 @@ class product_orthogonal(kernpart):
 
     def _set_params(self,x):
         """set the value of the parameters."""
-        self.k1._set_params(np.hstack((1.,x[1:self.k1.Nparam])))
-        self.k2._set_params(np.hstack((1.,x[self.k1.Nparam:])))
+        self.k1._set_params(x[:self.k1.Nparam])
+        self.k2._set_params(x[self.k1.Nparam:])
         self.params = x
 
     def _get_param_names(self):
         """return parameter names."""
-        return ['variance']+[self.k1.name + '_' + self.k1._get_param_names()[i+1] for i in range(self.k1.Nparam-1)] +  [self.k2.name + '_' + self.k2._get_param_names()[i+1] for i in range(self.k2.Nparam-1)]
-
+        return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k1._get_param_names()]
+    
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
         if X2 is None: X2 = X
@@ -71,7 +70,7 @@ class product_orthogonal(kernpart):
         target[0] += np.sum(K1*K2*partial)
         target[1:self.k1.Nparam] += self.params[0]* k1_target[1:]
         target[self.k1.Nparam:] += self.params[0]* k2_target[1:]
-        
+
     def dKdiag_dtheta(self,partial,X,target):
         """derivative of the diagonal of the covariance matrix with respect to the parameters."""
         target[0] += 1