0

For a homework problem we are asked to

Write a program that computes all the eigenvalues of the matrix A, using the Rayleigh Quotient iteration to find each eigenvalue.

i. Report the initial guess you used for each eigenvalue. Plot the error at the n iteration against the error at the n+1 iteration for each eigenvalue (in log log).

ii. Compute the Rayleigh quotient for x sampled on a discretization of the unit sphere (use spherical coordinates). Plot the result, for example with mesh (Θ,ϕ,r(Θ,ϕ)). Explain the number and location of the critical points.

I have completed the first part of the problem, but I am not sure I understand how to complete the second. After some research I have found various ways (Here and Here) to plot spherical coordinates in Python which has provided me some clues. Borrowing from those links and a few other sources I have

# The matrix for completness
A = np.matrix([[4,3,4], [3,1,3], [4,3,4]])
A = scipy.linalg.hessenberg(A)

num_points = 50
theta, phi = np.linspace(0, 2 * np.pi, num_points), np.linspace(0, np.pi, num_points)
x = np.sin(phi) * np.cos(theta)
y = np.sin(phi) * np.sin(theta)
z = np.cos(phi)

xyz = np.stack((x, y, z), axis=1)
rqs = np.array([np.dot(x, np.dot(A, np.conj(x).T)).item(0) for _,x in enumerate(xyz)])

In the first chunk of code (Θ,ϕ) are created and translated to Cartesian coordinates, effectively a discritization of the unit circle in the R^3. This allows me then to create the x vectors (xyz above) used to calculate the Rayleigh Quotients (rqs above) at each discritization point, the r(Θ,ϕ) of the spherical coordinates.

Though, now that I have the radial distance I am not sure how to again recreate x, y, z properly for a meshgrid to plot as a surface. This might be out of the realm of StackOverflow and more for Math.Stack, but I am also not sure if this plot should end up being a "warped plane" or a "warped sphere".

In this SO answer, linked above too, I think the answer lies though. In the chunk of code

theta, phi = np.linspace(0, 2 * np.pi, 40), np.linspace(0, np.pi, 40)
THETA, PHI = np.meshgrid(theta, phi)
R = np.cos(PHI**2)
X = R * np.sin(PHI) * np.cos(THETA)
Y = R * np.sin(PHI) * np.sin(THETA)
Z = R * np.cos(PHI)

R here I assume refers to the radial distance, though, this R is in a mesh when x, y, z are calculated. I have tried to reshape rqs above to the same dimension, but the values of the rqs do not line up with the subsequent grid and as a result produces obviously wrong plots.

I almost need a way to tie in the creation of the meshgrid with the calculation of the x. But the calculation seems to complex to be applied directly to the meshgrid..

How can I produce the spherical coordinates based plot given the radial distances?

EDIT: After some more searching I found this MatLab code which produces the desired plot, though I would still like to reproduce this in Python. I would like to say this MatLab code provides an overview of how to implement this in Python, but it seems to be some very old and esoteric code. Here is the plot(s) it produces

enter image description here

KDecker
  • 6,928
  • 8
  • 40
  • 81
  • What countour is it you want to plot? Is it the radius as a function of the angles? – JohanL Feb 26 '18 at 07:39
  • The contour is that of the Rayleigh Quotient of a Matrix `(M, x) = x* M x / x* x` (https://en.wikipedia.org/wiki/Rayleigh_quotient). It is a function of a matrix and a vector. Where the vector is a function of the angles; that is, the vectors are a discritization of the unit sphere in `R^3`. So the function should look like `r(M, x(Θ,ϕ)) = x(Θ,ϕ)* M x(Θ,ϕ) / x(Θ,ϕ)* x(Θ,ϕ)`. So then in the code above `x(Θ,ϕ) = xyz` and `r(M, x(Θ,ϕ)) = rqs`. The issue is I cannot figure out how to combine `rqs` with the creation of the angles and their subsequent translation to Cartesian coordinates. – KDecker Feb 26 '18 at 13:58
  • But `r` is a scalar value, that is a function of `Θ` and `ϕ`, while `M` can be considered constant (in this particular calculation)? – JohanL Feb 26 '18 at 14:57
  • Yes `r` will work out to be a scalar (The Rayleigh Quotient is a single scalar value). And I think this is where I am not sure, maybe need some Math.Stack in here. From how I understand it, I need to basically plot a unit sphere, as a surface plot, with a scalar value at each point on the sphere. So kind of a "warped sphere" so to say. `r` (well `rps`) represents the distance of each point on the unit sphere from the origin. – KDecker Feb 26 '18 at 15:08
  • Actually.. It should plot a unit sphere, but as a heat map. More intense colors at a given point on the unit sphere should represent larger values of `r`. I actually just found this MatLab code (https://people.eecs.berkeley.edu/~demmel/ma221_Fall04/Matlab/RayleighContour.m) which creates the plot needed. I will upload an image in the OP. – KDecker Feb 26 '18 at 15:16
  • Your Matlab code seems to be limited to show if it is larger or smaller than 1? – JohanL Feb 26 '18 at 15:36

1 Answers1

1

I don't know about the math or the Rayleigh Quotient but from what I gather, you want to calculate rqs as a function of the points of the unit sphere. For that, I would recommend to use meshgrid to generate value pairs for all Θ and ϕ. Then, since your matrix formula is defined for cartesian rather than spherical coordinates, I would convert my mesh to that and insert into the formula.

Then, finally, the result can be illustrated, using plot_surface on the unit sphere, where the (scaled) RQS data are used for facecolor:

import numpy as np
import scipy.linalg
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt

# The matrix for completness
A = np.matrix([[4,3,4], [3,1,3], [4,3,4]])
A = scipy.linalg.hessenberg(A)

num_points = 25
theta = np.linspace(0, 2 * np.pi, num_points)
phi = np.linspace(0, np.pi, num_points)

THETA, PHI = np.meshgrid(theta, phi)

X = np.sin(PHI) * np.cos(THETA)
Y = np.sin(PHI) * np.sin(THETA)
Z = np.cos(PHI)

# Calculate RQS for points on unit sphere:
RQS = np.empty(PHI.shape)
for theta_pos, phi_pos in itertools.product(range(num_points), range(num_points)):
    x = np.array([X[theta_pos, phi_pos],
                  Y[theta_pos, phi_pos],
                  Z[theta_pos, phi_pos]])
    RQS[theta_pos, phi_pos] = np.dot(x, np.dot(A, np.conj(x).T))

# normalize in range 0..1
maxRQS = abs(RQS).max()
N = (RQS+maxRQS)/(2*maxRQS)

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(
    X, Y, Z, rstride=1, cstride=1,
    facecolors=cm.jet(N),
    linewidth=0, antialiased=False, shade=False)
plt.show()

This gives the following result, which seems to agree with theMatlab contours plot in the OP.

enter image description here

Note that there may be better ways (vectorized) to convert the X, Y, and Z meshes to vectors, but the main execution time seems to be in the plottting anyway, so I won't spend the time trying to figure out exactly how.

JohanL
  • 6,671
  • 1
  • 12
  • 26