<filename>examples/_utilities/data_utilities.py
# SPDX-FileCopyrightText: Copyright 2021, <NAME> <<EMAIL>>
# SPDX-License-Identifier: BSD-3-Clause
# SPDX-FileType: SOURCE
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the license found in the LICENSE.txt file in the root
# directory of this source tree.


# =======
# Imports
# =======

import numpy


# ========================
# Generate Basis Functions
# ========================

def generate_basis_functions(n, m):
    """
    Generates a sample design matrix (also called basis functions).

    :param n: Nmber of rows of the design matrix.
    :type n: int

    :param m: Nmber of columns of the design matrix.
    :type m: int

    :return: Design matrix of the shape ``(n, m)``
    :rtype: numpy.ndarray

    The design matrix :math:`\\mathbf{X}` here is generated so that it is
    nearly singular. Using singular value decomposition, we generate
    :math:`\\mathbf{X}` by

    .. math::

        \\mathbf{X} = \\mathbf{U} \\mathbf{\\simga} \\mathbf{V}^{\\intercal},

    where :math:`\\mathbf{U}` and :math:`\\mathbf{V}` are orthonormal matrices
    and :math:`\\mathbf{\\sigma}` is the diagonal matrix of singular values.

    The orthonormal matrices are generated by Householer reflection

    .. math::

        \\mathbf{U} = \\mathbf{I} - 2 \\frac{ \\boldsymbol{u} \\otimes
        \\boldsymbol{u}}{\\| \\boldsymbol{u} \\|^2}

    and

    .. math::

        \\mathbf{V} = \\mathbf{I} - 2 \\frac{ \\boldsymbol{v} \\otimes
        \\boldsymbol{v}}{\\| \\boldsymbol{v} \\|^2}

    where :math:`\\boldsymbol{u}` and :math:`\\boldsymbol{v}` are random
    vectors. Also, the singular values (diagonals of :math:`\\mathbf{\\sigma}`)
    are produced by

    .. math::

        \\sigma_i = \\exp \\left( -40 \\left( \\frac{i}{m}
        \\right)^{\\frac{3}{4}} \\right).

    The reason to choose such function for singular values is that the
    generalized cross-validation function :math:`V(\\theta)`, which later we
    will plot, will have two non-trivial minima. The two minima is interesting
    to study for our application.
    """

    numpy.random.seed(31)

    if n > m:
        u = numpy.random.randn(n)
        U = numpy.eye(n) - 2.0 * numpy.outer(u, u) / numpy.linalg.norm(u)**2
        U = U[:, :m]
    else:
        u = numpy.random.randn(m)
        U = numpy.eye(m) - 2.0 * numpy.outer(u, u) / numpy.linalg.norm(u)**2
        U = U[:n, :]

    v = numpy.random.randn(m)
    V = numpy.eye(m) - 2.0 * numpy.outer(v, v.T) / numpy.linalg.norm(v)**2

    # sigma = numpy.exp(-20.0*(numpy.arange(m)/float(m))**(0.5))

    # good for n, m = 1000, 500
    sigma = numpy.exp(-40.0*(numpy.arange(m)/float(m))**(0.75))
    # sigma = numpy.exp(-20.0*(numpy.arange(m)/float(m))**(0.5)) * \
    #       numpy.sqrt(n/1000.0) * 1e2
    # sigma = numpy.exp(-10.0*(numpy.arange(m)/float(m))**(0.2)) * \
    #       numpy.sqrt(n/1000.0) * 1e2
    # sigma = numpy.sqrt(sigma)

    Sigma = numpy.diag(sigma)

    X = numpy.matmul(U, numpy.matmul(Sigma, V.T))

    return X


# ===================
# generate noisy data
# ===================

def generate_noisy_data(X, noise_level=4e-1):
    """
    Generates a column vector of noisy data by adding Gaussian noise to the
    mean function. The mean function is

    .. math::

        \\mu = \\mathbf{X} \\beta

    where :math:`\\mathbf{X}` is the design matrix and :math:`\\beta` is the
    coefficients of the linear function.

    :param X: Design matrix
    :type X: numpy.ndarray

    :param noise_level: Standard deviation of additive Gaussian noise.
    :type noise_level: float
    """

    # Size of basis function
    n, m = X.shape

    # beta = numpy.random.randn(m)
    beta = numpy.random.randn(m) / numpy.sqrt(n/1000.0)

    epsilon = noise_level * numpy.random.randn(n)

    # Data
    z = numpy.dot(X, beta) + epsilon

    return z


# ===============
# generate matrix
# ===============

def generate_matrix(n=1000, m=500, shift=1e-3):
    """
    Generates a symmetric positive semidefinite matrix.

    The matrix is generated by Grammian matrix as

    .. math::

        \\mathbf{K} =
        \\begin{cases}
            \\mathbf{X} \\mathbf{X}^{\\intercal} & n > m \\
            \\mathbf{X}^{\\intercal} \\mathbf{X} & n < m \\
        \\end{cases}

    where :math:`\\mathbf{X}` is of the size :math:`n \\times m`. Since
    :math:`\\mathbf{K}_0` is singular, we add a slight shift to its diagonals

    .. math::

        \\mathbf{K} = \\mathbf{K}_0 + s \\mathbf{I}.

    :param n: Number of rows of the design matrix.
    :type n: int

    :param m: Number of columns of the design matrix.
    :type m: int

    :param shift: A small number to be added to the diagonal elements of the
        matrix (in case the generated matrix is singular) to make the matrix
        slightly non-singular.
    :type shift: float
    """

    # Generate basis functions for linear regression
    X = generate_basis_functions(n, m)

    # Create Grammian matrix K0 from X
    if n > m:
        K0 = X.T.dot(X)

    else:
        K0 = X.dot(X.T)

    # Add shift to K0
    K = K0 + shift * numpy.eye(K0.shape[0], K0.shape[1])

    # Condition numbers
    cond_K0 = numpy.linalg.cond(K0)
    cond_K = numpy.linalg.cond(K)
    print('Cond K0: %0.2e, Cond K: %0.2e' % (cond_K0, cond_K))

    return K

<filename>examples/_utilities/data_utilities.py
# SPDX-FileCopyrightText: Copyright 2021, <NAME> <<EMAIL>>
# SPDX-License-Identifier: BSD-3-Clause
# SPDX-FileType: SOURCE
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the license found in the LICENSE.txt file in the root
# directory of this source tree.


# =======
# Imports
# =======

import numpy


# ========================
# Generate Basis Functions
# ========================

def generate_basis_functions(n, m):
    """
    Generates a sample design matrix (also called basis functions).

    :param n: Nmber of rows of the design matrix.
    :type n: int

    :param m: Nmber of columns of the design matrix.
    :type m: int

    :return: Design matrix of the shape ``(n, m)``
    :rtype: numpy.ndarray

    The design matrix :math:`\\mathbf{X}` here is generated so that it is
    nearly singular. Using singular value decomposition, we generate
    :math:`\\mathbf{X}` by

    .. math::

        \\mathbf{X} = \\mathbf{U} \\mathbf{\\simga} \\mathbf{V}^{\\intercal},

    where :math:`\\mathbf{U}` and :math:`\\mathbf{V}` are orthonormal matrices
    and :math:`\\mathbf{\\sigma}` is the diagonal matrix of singular values.

    The orthonormal matrices are generated by Householer reflection

    .. math::

        \\mathbf{U} = \\mathbf{I} - 2 \\frac{ \\boldsymbol{u} \\otimes
        \\boldsymbol{u}}{\\| \\boldsymbol{u} \\|^2}

    and

    .. math::

        \\mathbf{V} = \\mathbf{I} - 2 \\frac{ \\boldsymbol{v} \\otimes
        \\boldsymbol{v}}{\\| \\boldsymbol{v} \\|^2}

    where :math:`\\boldsymbol{u}` and :math:`\\boldsymbol{v}` are random
    vectors. Also, the singular values (diagonals of :math:`\\mathbf{\\sigma}`)
    are produced by

    .. math::

        \\sigma_i = \\exp \\left( -40 \\left( \\frac{i}{m}
        \\right)^{\\frac{3}{4}} \\right).

    The reason to choose such function for singular values is that the
    generalized cross-validation function :math:`V(\\theta)`, which later we
    will plot, will have two non-trivial minima. The two minima is interesting
    to study for our application.
    """

    numpy.random.seed(31)

    if n > m:
        u = numpy.random.randn(n)
        U = numpy.eye(n) - 2.0 * numpy.outer(u, u) / numpy.linalg.norm(u)**2
        U = U[:, :m]
    else:
        u = numpy.random.randn(m)
        U = numpy.eye(m) - 2.0 * numpy.outer(u, u) / numpy.linalg.norm(u)**2
        U = U[:n, :]

    v = numpy.random.randn(m)
    V = numpy.eye(m) - 2.0 * numpy.outer(v, v.T) / numpy.linalg.norm(v)**2

    # sigma = numpy.exp(-20.0*(numpy.arange(m)/float(m))**(0.5))

    # good for n, m = 1000, 500
    sigma = numpy.exp(-40.0*(numpy.arange(m)/float(m))**(0.75))
    # sigma = numpy.exp(-20.0*(numpy.arange(m)/float(m))**(0.5)) * \
    #       numpy.sqrt(n/1000.0) * 1e2
    # sigma = numpy.exp(-10.0*(numpy.arange(m)/float(m))**(0.2)) * \
    #       numpy.sqrt(n/1000.0) * 1e2
    # sigma = numpy.sqrt(sigma)

    Sigma = numpy.diag(sigma)

    X = numpy.matmul(U, numpy.matmul(Sigma, V.T))

    return X


# ===================
# generate noisy data
# ===================

def generate_noisy_data(X, noise_level=4e-1):
    """
    Generates a column vector of noisy data by adding Gaussian noise to the
    mean function. The mean function is

    .. math::

        \\mu = \\mathbf{X} \\beta

    where :math:`\\mathbf{X}` is the design matrix and :math:`\\beta` is the
    coefficients of the linear function.

    :param X: Design matrix
    :type X: numpy.ndarray

    :param noise_level: Standard deviation of additive Gaussian noise.
    :type noise_level: float
    """

    # Size of basis function
    n, m = X.shape

    # beta = numpy.random.randn(m)
    beta = numpy.random.randn(m) / numpy.sqrt(n/1000.0)

    epsilon = noise_level * numpy.random.randn(n)

    # Data
    z = numpy.dot(X, beta) + epsilon

    return z


# ===============
# generate matrix
# ===============

def generate_matrix(n=1000, m=500, shift=1e-3):
    """
    Generates a symmetric positive semidefinite matrix.

    The matrix is generated by Grammian matrix as

    .. math::

        \\mathbf{K} =
        \\begin{cases}
            \\mathbf{X} \\mathbf{X}^{\\intercal} & n > m \\
            \\mathbf{X}^{\\intercal} \\mathbf{X} & n < m \\
        \\end{cases}

    where :math:`\\mathbf{X}` is of the size :math:`n \\times m`. Since
    :math:`\\mathbf{K}_0` is singular, we add a slight shift to its diagonals

    .. math::

        \\mathbf{K} = \\mathbf{K}_0 + s \\mathbf{I}.

    :param n: Number of rows of the design matrix.
    :type n: int

    :param m: Number of columns of the design matrix.
    :type m: int

    :param shift: A small number to be added to the diagonal elements of the
        matrix (in case the generated matrix is singular) to make the matrix
        slightly non-singular.
    :type shift: float
    """

    # Generate basis functions for linear regression
    X = generate_basis_functions(n, m)

    # Create Grammian matrix K0 from X
    if n > m:
        K0 = X.T.dot(X)

    else:
        K0 = X.dot(X.T)

    # Add shift to K0
    K = K0 + shift * numpy.eye(K0.shape[0], K0.shape[1])

    # Condition numbers
    cond_K0 = numpy.linalg.cond(K0)
    cond_K = numpy.linalg.cond(K)
    print('Cond K0: %0.2e, Cond K: %0.2e' % (cond_K0, cond_K))

    return K

# SPDX-FileCopyrightText: Copyright 2021, <NAME> <<EMAIL>> # SPDX-License-Identifier: BSD-3-Clause # SPDX-FileType: SOURCE # # This program is free software: you can redistribute it and/or modify it # under the terms of the license found in the LICENSE.txt file in the root # directory of this source tree. # ======= # Imports # ======= # ======================== # Generate Basis Functions # ======================== Generates a sample design matrix (also called basis functions).      :param n: Nmber of rows of the design matrix.     :type n: int      :param m: Nmber of columns of the design matrix.     :type m: int      :return: Design matrix of the shape ``(n, m)``     :rtype: numpy.ndarray      The design matrix :math:`\\mathbf{X}` here is generated so that it is     nearly singular. Using singular value decomposition, we generate     :math:`\\mathbf{X}` by      .. math::          \\mathbf{X} = \\mathbf{U} \\mathbf{\\simga} \\mathbf{V}^{\\intercal},      where :math:`\\mathbf{U}` and :math:`\\mathbf{V}` are orthonormal matrices     and :math:`\\mathbf{\\sigma}` is the diagonal matrix of singular values.      The orthonormal matrices are generated by Householer reflection      .. math::          \\mathbf{U} = \\mathbf{I} - 2 \\frac{ \\boldsymbol{u} \\otimes         \\boldsymbol{u}}{\\| \\boldsymbol{u} \\|^2}      and      .. math::          \\mathbf{V} = \\mathbf{I} - 2 \\frac{ \\boldsymbol{v} \\otimes         \\boldsymbol{v}}{\\| \\boldsymbol{v} \\|^2}      where :math:`\\boldsymbol{u}` and :math:`\\boldsymbol{v}` are random     vectors. Also, the singular values (diagonals of :math:`\\mathbf{\\sigma}`)     are produced by      .. math::          \\sigma_i = \\exp \\left( -40 \\left( \\frac{i}{m}         \\right)^{\\frac{3}{4}} \\right).      The reason to choose such function for singular values is that the     generalized cross-validation function :math:`V(\\theta)`, which later we     will plot, will have two non-trivial minima. The two minima is interesting     to study for our application. # sigma = numpy.exp(-20.0*(numpy.arange(m)/float(m))**(0.5)) # good for n, m = 1000, 500 # sigma = numpy.exp(-20.0*(numpy.arange(m)/float(m))**(0.5)) * \ #       numpy.sqrt(n/1000.0) * 1e2 # sigma = numpy.exp(-10.0*(numpy.arange(m)/float(m))**(0.2)) * \ #       numpy.sqrt(n/1000.0) * 1e2 # sigma = numpy.sqrt(sigma) # =================== # generate noisy data # =================== Generates a column vector of noisy data by adding Gaussian noise to the     mean function. The mean function is      .. math::          \\mu = \\mathbf{X} \\beta      where :math:`\\mathbf{X}` is the design matrix and :math:`\\beta` is the     coefficients of the linear function.      :param X: Design matrix     :type X: numpy.ndarray      :param noise_level: Standard deviation of additive Gaussian noise.     :type noise_level: float # Size of basis function # beta = numpy.random.randn(m) # Data # =============== # generate matrix # =============== Generates a symmetric positive semidefinite matrix.      The matrix is generated by Grammian matrix as      .. math::          \\mathbf{K} =         \\begin{cases}             \\mathbf{X} \\mathbf{X}^{\\intercal} & n > m \\             \\mathbf{X}^{\\intercal} \\mathbf{X} & n < m \\         \\end{cases}      where :math:`\\mathbf{X}` is of the size :math:`n \\times m`. Since     :math:`\\mathbf{K}_0` is singular, we add a slight shift to its diagonals      .. math::          \\mathbf{K} = \\mathbf{K}_0 + s \\mathbf{I}.      :param n: Number of rows of the design matrix.     :type n: int      :param m: Number of columns of the design matrix.     :type m: int      :param shift: A small number to be added to the diagonal elements of the         matrix (in case the generated matrix is singular) to make the matrix         slightly non-singular.     :type shift: float # Generate basis functions for linear regression # Create Grammian matrix K0 from X # Add shift to K0 # Condition numbers