FreeConvolutionCode/test_multiplicative_convolution.py at main · Alice94/FreeConvolutionCode · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
import spectral_measure as sm
import numpy
import matplotlib.pyplot as plt
import freeconvolution
import scipy
from random import choices

# Set up parameters
do_plots = 1
epsilon = 0.05
matrix_size = 5000
bin_number = int(100)

for n_example in range(0, 5):

  mu1 = sm.spectral_measure()
  mu2 = sm.spectral_measure()

  if (n_example == 0):
    param = 0.1
    # ~ param2 = 0.2
    a2 = 1
    b2 = 3
    mu1.set_marchenko_pastur(param)
    mu2.set_uniform(a2, b2)
    N1 = 400
    N2 = 400
    N = 4000
    m = 20
    name = "MPxunif.eps"
    A1 = mu1.get_random_matrix(matrix_size)
    A2 = mu2.get_random_matrix(matrix_size)
    A = numpy.matmul(A1, A2)
    eigenvalues = scipy.linalg.eigvals(A)
    print(numpy.min(eigenvalues), numpy.max(eigenvalues))

  elif (n_example == 1):
    N1 = 4000
    N2 = 4000
    N = 4000
    m = 20
    param = 0.2
    mu1.set_shifted_semicircle(1)
    mu2.set_marchenko_pastur(param)
    name = "SemicirclexMP0.2.eps"
    A1 = mu1.get_random_matrix(matrix_size)
    A2 = mu2.get_random_matrix(matrix_size)
    A = numpy.matmul(A1, A2)
    eigenvalues = scipy.linalg.eigvals(A)
    epsilon = 0.05

  elif (n_example == 2):
    N1 = 400
    N2 = 400
    N = 400
    m = 20
    param = 1
    mu1.set_shifted_semicircle(1)
    mu2.set_shifted_semicircle(1)
    name = "Semicirclexsemicircle.eps"
    A1 = mu1.get_random_matrix(matrix_size)
    A2 = mu2.get_random_matrix(matrix_size)
    A = numpy.matmul(A1, A2)
    eigenvalues = scipy.linalg.eigvals(A)

  elif (n_example == 3):
    N1 = 400
    N2 = 4000
    N = 4000
    m = 24
    mu1.set_shifted_semicircle(1)
    weights = numpy.array(numpy.ones(7))
    weights = weights / sum(weights)
    points = numpy.array(numpy.linspace(1,4,num=7))
    mu2.set_discrete(points, weights)
    name = "semicirclexdiscrete.eps"
    A1 = mu1.get_random_matrix(matrix_size)
    A2 = mu2.get_random_matrix(matrix_size)
    A = numpy.matmul(A1, A2)
    eigenvalues = scipy.linalg.eigvals(A)
    epsilon = 0.08

  elif (n_example == 4):
    a1 = 1
    a2 = 1
    b1 = 2
    b2 = 3

    mu1.set_uniform(a1, b1)
    mu2.set_uniform(a2, b2)
    N1 = 4000
    N2 = 4000
    N = 4000
    m = 20
    name = "Unif1xUnif2.eps"
    A1 = mu1.get_random_matrix(matrix_size)
    A2 = mu2.get_random_matrix(matrix_size)
    A = numpy.matmul(A1, A2)
    eigenvalues = scipy.linalg.eigvals(A)

  [a_prod, b_prod, t, approx_mu, mutimes] = freeconvolution.free_multiplicative_convolution(mu1, mu2, N1, N2, N, m, epsilon, do_plots)
  plt.subplot(4, 3, 10)
  (counts, bins) = numpy.histogram(eigenvalues, bins=bin_number)
  plt.hist(bins[:-1], bins, weights=bin_number / (matrix_size * (mutimes.b - mutimes.a)) * counts)
  plt.plot(t, approx_mu, 'r')

  plt.subplot(4, 3, 11)
  plt.plot(t, abs(approx_mu), 'r')

  figure = plt.gcf() # get current figure
  figure.set_size_inches(16, 12)
  plt.savefig(name, dpi=10000)
  # ~ plt.show()

  del(mu1)
  del(mu2)