Hi Sebastien,

thanks for your reply! I admit, in the scattering context the assumption of "lossless" transmission doesn't make sense. Parametrization to a fixed transmission is surely a better idea, especially considering that this can be confirmed/measured in an experimental setup.

Cheers,
Matthias

Hi Sebastien,
thanks a lot for your reply! You're right, the underlying problem was my physical mis-interpretation of the scattering matrix. And with the correction you stated it of course works out to the quarter circle law in this specific case, independent of the variance of distribution.

One more thought on the normalization of the transmission matrix:
Given the situation, when we illuminate all N input pixels with unity amplitude, implement a given input mode n as a phase function and measure the response at all M output modes - the total recorded energy will be equal to the input (=N*1), or less. So, the mean of the M output pixel amplitudes will be <= N/M.

So, this measurement basically reflects one column of the transmission matrix (that is, if we measure amplitude and phase). And it would be the same for all other column measurements - the mean absolute value should be <= N/M = gamma.

So, could it be, that for the case of zero absorption we have to scale to mean(abs(k_ij)) = gamma? And for practical experiments the elements will of course have smaller absolute values.

My updated code looks as follows:

#%% calculations
from numpy import sqrt, angle, mean, sum
from numpy.random import randn
from scipy.linalg import svd
import matplotlib.pyplot as plt
m = 1024
gamma = 2
n = gamma * m

K = randn(m,n)+ 1j*randn(m,n) # transmission matrix
normFactor = mean(mean(abs(K))) # energy conservation (here zero absorption)
K = K/normFactor*gamma # normalization of transmission matrix
U,s,V = svd(K)

# singular value normalization
sNorm = s*sqrt(m)/sqrt(sum(s**2)) # as in Popoff et al. --> corrected

#%% plotting
plt.figure()
plt.hist(sNorm, bins = 50, normed = True)
plt.title('histogram of singular values, normalization sqrt(m)/sqrt(sum(s**2))')
plt.xlabel('normalized singular value')
plt.ylabel('probability density')

plt.figure()
plt.subplot(1,2,1)
plt.imshow(abs(K),cmap='gray')
plt.colorbar(fraction = 0.046, pad=0.04)
plt.title('Transmission matrix in amplitude and phase')

plt.subplot(1,2,2)
plt.imshow(angle(K),cmap='jet')
plt.colorbar(fraction = 0.046, pad=0.04)

plt.show()

In this case mean(mean(abs(K))) yields 2=gamma, which I would expect assuming "zero loss".

Cheers,
Matthias

Hi community,

I'm working on incorporating a model of a scattering medium into a simulation and I need to check, whether my model is correct. It is assumed that the conditions for random matrix theory to be a valid description are met. I apply SVD to compare to the results of Popoff et al. (doi.org/10.1088/1367-2630/13/12/123021). For SVD implementation I relied on the paper by Edelman&Wang (web.eecs.umich.edu/~rajnrao/Acta05rmt.pdf). My problem lies in the normaliziation of the singular values.

I assume normal distribution N(0,1) for both real and imaginary part of the transmission matrix. The matrix is normalized to the biggest absolute value, as I assume a passive transmission medium with no amplification possible.
Now, the quarter circle law distribution for the singular values is to be expected for this case of gamma=1. The probability distribution follows this curve, but I obtain an incorrect range for the singular values (quarter circle law predicts [0,2]).
After the scaling of K, the variance its elements will not be 1 anymore. So, I found the previously applied scaling factor, divided by sqrt(2), to be a working conversion to the [0,2] interval predicted by the quarter circle law.

However, the physical interpretation is not clear to me. By visualizing the singular value I want to verify, that I use a description previously shown to be valid for certain scattering media.

Would someone know, how to correctly scale the singular values in order to show, that the example follows the quarter circle law?
I'd be grateful for any hint.
Cheers, Matthias

Here's the numeric Python code I used:

from numpy import sqrt, angle, amax, sum
from numpy.random import randn
from scipy.linalg import svd
import matplotlib.pyplot as plt
m = 1024
gamma = 1
n = gamma * m

K = randn(m,n)+ 1j*randn(m,n) # transmission matrix
normFactor = amax(amax(abs(K))) # biggest amplitude element in K
K = K/normFactor # normalization
U,s,V = svd(K)

# different normalizations
sNorm1 = s/sqrt(m) # as in Edelman & Wang
sNorm2 = s/sqrt(sum(s**2)) # as in Popoff et al.
sNorm3 = normFactor/sqrt(2)*s/(sqrt(m)) # results in [0,2] interval (as expected)

#%% plotting
plt.figure()
plt.hist(sNorm1, bins = 50, normed = True)
plt.title('histogram of singular values, normalization 1/sqrt(m)')
plt.xlabel('normalized singular value')
plt.ylabel('probability density')

plt.figure()
plt.hist(sNorm2, bins = 50, normed = True)
plt.title('histogram of singular values, normalization 1/sqrt(sum(s**2))')
plt.xlabel('normalized singular value')
plt.ylabel('probability density')

plt.figure()
plt.hist(sNorm3, bins = 50, normed = True)
plt.title('histogram of singular values, normalization normFactor/sqrt(2)*/(sqrt(m))')
plt.xlabel('normalized singular value')
plt.ylabel('probability density')

plt.show()