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.
getSingularValues()	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()