Orthogonal initialization of weight matrix

Question

Searching for the way to initialize the matrix weights as orthogonal (i.e. W*W^T = I and all the eigenvalues are equal either 1 or -1),(I was wrong) I found this post with a code from Lasagne:

https://stats.stackexchange.com/questions/228704

However, this code does not produce orthogonal, but unitary matrix with complex eigenvalues, which still satisfies equality W*W^T = I quite well:

import numpy as np


np.random.seed(1000)


def floatX(arr):
    """Converts data to a numpy array of dtype ``theano.config.floatX``.
    Parameters
    ----------
    arr : array_like
        The data to be converted.
    Returns
    -------
    numpy ndarray
        The input array in the ``floatX`` dtype configured for Theano.
        If `arr` is an ndarray of correct dtype, it is returned as is.
    """
    return np.asarray(arr, dtype=np.float32)

def sample(shape,gain=1.0):

        flat_shape = (shape[0], np.prod(shape[1:]))
        a = np.random.normal(0, 1, flat_shape).astype(np.float32)  #get_rng().normal(0.0, 1.0, flat_shape)
        u, _, v = np.linalg.svd(a, full_matrices=False)
        # pick the one with the correct shape
        q = u if u.shape == flat_shape else v
        q = q.reshape(shape)
        res = floatX(gain * q)

        #eigenvalues:
        eigenvalues, _ = np.linalg.eig(u)

        print(str(np.min(eigenvalues))+" "+ str(eigenvalues[1024])+ " "+ str(np.max(eigenvalues))  )

        return res 


W = sample([2048, 2048])


print(np.linalg.norm(W.dot(np.transpose(W))-np.eye(2048)))

The output is

(-1+0j) (0.376566+0.92639j) (1+0j)

for some of the eigenvalues and

1.00117376019e-05

for the proximity of the W*W^T and I.

1) Is this a bug in Lasagne, as numpy.linalg.svd indeed produces unitary matrix as it is stated in its manual, and

2) How can I initialize a (non-square) orthogonal matrix for myself in Python, except the unit one?

Thank you!

EDIT: Thanks everyone! From all the explanations I've managed to create a simple Python solution:

import numpy as np

## Using a subset of vectors
np.random.seed(1000)

p = 4096
A = np.random.normal(0,1,(p,p))

q, _ = np.linalg.qr(A)  # q -orthonormal

print(q.shape) #(4096,4096)

q = q[:,::2] #cut q

print(q.shape) #(4096,2048) m>n ==> q^T * q = I_n

print(np.linalg.norm(np.transpose(q).dot(q)-np.eye(2048))) # 4.26e-14

Answer 1

How can I initialize a (non-square) orthogonal matrix for myself in Python, except the unit one?

One way to make the orthogonal matrix is to get the subset of eigenvectors of any positive definite matrix.

Here's Python code:

import numpy as np
import math
from scipy.linalg import toeplitz
from numpy import linalg as LA

# create the positive definite matrix P
p = 4
A = np.random.rand(p,p)
P = (A + np.transpose(A))/2 + p*np.eye(p)

# get the subset of its eigenvectors
vals, vecs = LA.eig(P)
w = vecs[:,0:p-1]

#check that it's orthogonal 
print("matrix shape: ",w.shape)
print("check orthogonality: ",np.matmul(np.transpose(w),w))

matrix shape:  (4, 3)
check orthogonality:  [[  1.00000000e+00   3.88578059e-16  -6.21031004e-16]
 [  3.88578059e-16   1.00000000e+00   5.55111512e-16]
 [ -6.21031004e-16   5.55111512e-16   1.00000000e+00]]

UPDATE

This matrix does not (and cannot) have complex elements, of course. Here's the proof:

Python Code:

import numpy as np
import math
from scipy.linalg import toeplitz
from numpy import linalg as LA

# create the positive definite matrix P
p = 400
A = np.random.rand(p,p)
P = (A + np.transpose(A))/2 + p*np.eye(p)

# get the subset of its eigenvectors
vals, vecs = LA.eig(P)
w = vecs[:,0:p]

#check that it's orthogonal 
print("matrix shape: ",w.shape)
print("check orthogonality: ",np.matmul(np.transpose(w),w))
print("has complex elements: ",np.max(np.iscomplex(w)))

Output:

matrix shape:  (400, 400)
check orthogonality:  [[  1.00000000e+00   1.33573708e-16   1.67509236e-16 ...,  -5.20417043e-18
    7.50267903e-17  -8.84708973e-17]
 [  1.33573708e-16   1.00000000e+00  -7.39425882e-15 ...,  -1.08940634e-15
    1.03172679e-15  -6.50521303e-16]
 [  1.67509236e-16  -7.39425882e-15   1.00000000e+00 ...,   3.15318518e-15
    1.58553726e-15   1.27155231e-15]
 ..., 
 [ -5.20417043e-18  -1.08940634e-15   3.15318518e-15 ...,   1.00000000e+00
    4.93181884e-15  -1.13573213e-13]
 [  7.50267903e-17   1.03172679e-15   1.58553726e-15 ...,   4.93181884e-15
    1.00000000e+00  -9.91155075e-13]
 [ -8.84708973e-17  -6.50521303e-16   1.27155231e-15 ...,  -1.13573213e-13
   -9.91155075e-13   1.00000000e+00]]
has complex elements:  False