blindem_linear_operators.day_by_day.py

import random

import numpy
import numpy.matlib

N = 10
S = 3

Omega = random.sample(range(N), S)
s = numpy.matlib.randn(S,1)

C_Omegat = numpy.zeros([N,S])
C_Omegat[Omega,:] = numpy.eye(S)

x = C_Omegat@s


### Day 2

import numpy as np

L = 4

def make_connected_erdos_rendi(N=10, p=0.1):
  while True:
    S = np.array(np.random.rand(N, N) < p, dtype=float)
    S = np.triu(S, k=1)
    S = S + S.T
    # Make sure it is connected.
    L = np.diag(np.sum(S, 0)) - S
    eigvals = np.linalg.eig(L)[0]
    if sum(abs(eigvals) < 1e-6) == 1:
      return S

h = np.random.randn(L, 1)
x = np.random.randn(N, 1)

S = make_connected_erdos_rendi()

eigvals, V = np.linalg.eig(S)
V = np.matrix(V)
U = np.linalg.inv(V)
Psi = np.vander(eigvals, len(h), increasing=True)
Psi = np.matrix(Psi)

#print((U@C_Omegat).shape)

# project
Omegaperp = np.setdiff1d(range(N), Omega, assume_unique=False)

import matplotlib.pyplot

Uprojected = U
Uprojected[:,Omegaperp] = 0
matplotlib.pyplot.spy(Uprojected)
print(Uprojected)


### Day 3

def linop_M_Omega(Z):
    ret = np.zeros(N, dtype=np.csingle)
    for n in range(N):
        #print(U_Omega[n,:].shape, Z.shape, np.atleast_2d(Psi[n,:]).T.shape)
        ret[n] = U_Omega[n,:]@Z@np.atleast_2d(Psi[n,:]).T
    return ret

a = linop_M_Omega(np.random.rand(N,L))
print(np.abs(a))

### Day 4

def linop_M_Omega_adj(z):
    R = np.zeros([N, L], dtype=np.csingle)
    for n in range(N):
        R += z[n]*(np.conjugate(np.atleast_2d(U_Omega[n,:]).T)@np.conjugate(Psi[n,:]))
    return R

#B = linop_M_Omega_adj(b)

#(M^*_{i \Omega_i} M_{i \Omega_i} - I d
#print(linop_M_Omega_adj(linop_M_Omega(np.random.rand(N,L))) - np.eye(N=N,M=L))
#matplotlib.pyplot.spy(linop_M_Omega_adj(linop_M_Omega(np.random.rand(N,L))))
#matplotlib.pyplot.spy(linop_M_Omega_adj(linop_M_Omega(np.random.rand(N,L))) - np.eye(N=N,M=L))

# TODO Consider randomize to form A?
A = U@C_Omegat

def linop_A(Z):
    r = np.zeros(N, dtype=np.csingle)
    for n in range(N):
        r[n] = A[n,:]*Z*np.atleast_2d(Psi[n,:]).T
    return r

def linop_A_adj(z):
    R = np.zeros([S, L], dtype=np.csingle)
    for n in range(N):
        R += z[n]*(np.conjugate(np.atleast_2d(A[n,:]).T)@np.conjugate(Psi[n,:]))
    return R

# FIXME linop_M_Omega and linop_A are very similar, only the matrix U or A change, share code parameterizing.

blindem_linear_operators.py

import random

import numpy
import numpy.matlib

L = 4
N = 10
S = 3

# Sample S unique elements from 1,...,N.
Omega = random.sample(range(N), S)
s = numpy.matlib.randn(S,1)

C_Omegat = numpy.zeros([N,S])
C_Omegat[Omega,:] = numpy.eye(S)

x = C_Omegat@s

import numpy as np

def make_connected_erdos_rendi(N=10, p=0.1):
  while True:
    S = np.array(np.random.rand(N, N) < p, dtype=float)
    S = np.triu(S, k=1)
    S = S + S.T
    # Make sure it is connected.
    L = np.diag(np.sum(S, 0)) - S
    eigvals = np.linalg.eig(L)[0]
    if sum(abs(eigvals) < 1e-6) == 1:
      return S

#h = np.random.randn(L, 1)
x = np.random.randn(N, 1)

GSO = make_connected_erdos_rendi(N=N)

eigvals, V = np.linalg.eig(GSO)
V = np.matrix(V)
U = np.linalg.inv(V)
Psi = np.vander(eigvals, L, increasing=True)
Psi = np.matrix(Psi)

# project
Omegaperp = np.setdiff1d(range(N), Omega, assume_unique=False)

import matplotlib.pyplot

U_Omega = U
U_Omega[:,Omegaperp] = 0 # ojo!
#matplotlib.pyplot.spy(U_Omega);

def linop_M_Omega(Z):
    r = np.zeros(N, dtype=np.csingle)
    for n in range(N):
        #print(U_Omega[n,:].shape, Z.shape, np.atleast_2d(Psi[n,:]).T.shape)
        r[n] = (U_Omega[n,:]@Z)@np.atleast_2d(Psi[n,:]).T
    return r

#b = linop_M_Omega(np.random.rand(N,L))

def linop_M_Omega_adj(z):
    R = np.zeros([N, L], dtype=np.csingle)
    for n in range(N):
        R += z[n]*(np.conjugate(np.atleast_2d(U_Omega[n,:]).T)@np.conjugate(Psi[n,:]))
    return R

#B = linop_M_Omega_adj(b)

A = U@C_Omegat

def linop_A(Z):
    r = np.zeros(N, dtype=np.csingle)
    for n in range(N):
        r[n] = A[n,:]@Z@np.atleast_2d(Psi[n,:]).T
    return r

def linop_A_adj(z):
    R = np.zeros([S, L], dtype=np.csingle)
    for n in range(N):
        R += z[n]*(np.conjugate(np.atleast_2d(A[n,:]).T)@np.conjugate(Psi[n,:]))
    return R

#(M^*_{i \Omega_i} M_{i \Omega_i} - Id
#(A^*_i A - Id)

Z = np.random.rand(N,L)
X = linop_M_Omega_adj(linop_M_Omega(Z))
Y = linop_A_adj(linop_A(C_Omegat.T@Z))
n1 = np.linalg.norm(X - np.eye(N=N,M=L), ord=2)
n2 = np.linalg.norm(Y - np.eye(N=S,M=L), ord=2)

print('Support=', Omega, '\n')

print('X=\n',X,'\n')
print('Y=\n',Y,'\n')

print('X-I=\n',X - np.block([C_Omegat, np.zeros([N,L-S])]),'\n')
print('Y-I=\n',Y - np.eye(N=S,M=L),'\n')

print(np.linalg.norm(np.block([C_Omegat, np.zeros([N,L-S])]), ord=2))
print(np.linalg.norm(np.eye(N=S,M=L), ord=2))

x is sparse with N elements and S nnz

x = C_Omegat s
[N] = [NxS] [S]