The Singular Value Decomposition (SVD)¶
PyA’s SVD class provides the means to carry out a singular value decomposition as is, e.g., useful to recover line profiles.
- class PyAstronomy.pyasl.SVD¶
Use Singular Value Decomposition (SVD) to obtain broadening function.
The technique is, e.g., described in Rucinski 1999, “Determination of Broadening Functions Using the Singular-Value Decomposition (SVD) Technique” and references therein.
Methods
decompose
(template, bn)Carry out the SVD of the "design matrix".
getBroadeningFunction
(flux[, wlimit, asarray])Calculate the broadening function.
getModel
(broad[, asarray, modelIndices])Calculates a model resulting from template and broadening function.
getRVAxis
(binsize, refWvl)Calculate a radial velocity axis for the broadening function.
Access the singular values.
- decompose(template, bn)¶
Carry out the SVD of the “design matrix”.
This method creates the “design matrix” by applying a bin-wise shift to the template and uses numpy’s svd algorithm to carry out the decomposition.
The design matrix, des is written in the form: “des = u * w * transpose(v)”. The matrices des, w, and u are stored in homonymous attributes.
- Parameters
- templatearray or matrix
The template with respect to which the broadening function shall be calculated.
- bnint
Width (number of elements) of the broadening function. Needs to be odd.
- getBroadeningFunction(flux, wlimit=0.0, asarray=False)¶
Calculate the broadening function.
Note
The decompose method has to be called first. On this call, the template is specified.
- Parameters
- fluxarray or matrix
The observed function (flux).
- wlimitfloat, optional
A limit to be applied to the singular values. Values smaller than the specified limit will be discarded in calculating the broadening function.
- asarraybool, optional
If True, the broadening function will be returned as an array instead of a matrix. The default is False.
- Returns
- Broadening functionmatrix, array
The broadening function. The return type can be set by the asarray flag.
- getModel(broad, asarray=False, modelIndices=False)¶
Calculates a model resulting from template and broadening function.
- Parameters
- broadarray or matrix
The broadening function
- asarraybool, optional
If True, the broadening function will be returned as an array instead of a matrix. The default is False.
- modelIndicesbool, optional
If True, the method returns also an array of indices, which refer to the template, for which the model is calculated. Note that the model cannot be calculated at the edges. The default is False.
- getRVAxis(binsize, refWvl)¶
Calculate a radial velocity axis for the broadening function.
Note
Here, it is assumed that the broadening function refers to a spectrum in wavelength space.
- Parameters
- binsizefloat
Size of spectral bins.
- refWvlfloat
The reference wavelength.
- Returns
- Radial velocity axisarray
An array containing the radial velocity shift pertaining to the individual bins in the broadening function. The unit is km/s.
- getSingularValues()¶
Access the singular values.
- Returns
- Singular valuesmatrix
The singular values.
Example: Delta functions and rotational broadening¶
import numpy as np
import matplotlib.pylab as plt
from PyAstronomy import pyasl
# Get some "data"
wvl = np.arange(5000., 5010., 0.02)
template = np.ones(len(wvl))
# There are two sharp lines in the template
template[100] = 0.5
template[115] = 0.3
# Apply rotational broadening to the delta spectrum
nflux = pyasl.rotBroad(wvl, template, 0.5, 23.45, edgeHandling="firstlast")
# Carry out decomposition
svd = pyasl.SVD()
# Use 51 bins for the width of the broadening function.
# Needs to be an odd number.
svd.decompose(template, 51)
# Obtain the broadening function needed to
# recover "observed" spectrum. Note that the
# edges (51/2 bins) will actually be neglected.
b = svd.getBroadeningFunction(nflux)
# Get the model, which results from the broadening
# function and the template; obtain the indices
# where it applies, too.
m, mind = svd.getModel(b, modelIndices=True)
# Plot the outcome
plt.plot(b, 'bp-')
plt.plot(mind, m, 'r.')
plt.plot(nflux, 'g--')
plt.show()
Example: Adding noise and neglecting singular values¶
import numpy as np
import matplotlib.pylab as plt
from PyAstronomy import pyasl
# Get some "data"
wvl = np.arange(5000., 5010., 0.02)
template = np.ones(len(wvl))
# There are two sharp lines in the template
template[100] = 0.5
template[115] = 0.3
# Apply rotational broadening to the delta spectrum
nflux = pyasl.rotBroad(wvl, template, 0.5, 23.45, edgeHandling="firstlast")
nflux += np.random.normal(0., 0.005, len(nflux))
# Carry out decomposition
svd = pyasl.SVD()
svd.decompose(template, 51)
# Access the singular values
sv = svd.getSingularValues()
# Calculate the reduced chi square as a function of the number
# of singular values neglected in the calculation of the
# model.
chi = []
for i in range(1, len(sv), 5):
b = svd.getBroadeningFunction(nflux, wlimit=sorted(sv)[i])
m, mind = svd.getModel(b, modelIndices=True, asarray=True)
chi.append(((nflux[mind] - m)**2/0.005**2).sum() / len(mind))
plt.title("Reduced $\chi^2$ vs. number of neglected singular values")
plt.plot(range(1, len(sv), 5), chi, 'bp-')
plt.show()