text:
BayesOpt
Bayesian optimization

# Introduction to Bayesian Optimization

Many problems in engineering, computer science, economics, etc., require to find the extremum of a real valued function. These functions are typically continuous and sometimes smooth (e.g.: Lipschitz continuous). However, those functions do not have a closed-form expression or might be multimodal, where some of the local extrema might have a bad outcome compared to the global extremum. The evaluation of those functions might be costly.

Global optimization is a special case of non-convex optimization where we want to find the global extremum of a real valued function, that is, the target function. The search is done by some pointwise evaluation of the target function.

The objective of a global optimization algorithm is to find the sequence of points

which converges to the point , that is, the extremum of the target function, when is large. The algorithm should be able to find that sequence at least for all functions from a given family.

As explained in[18], this search procedure is a sequential decision making problem where point at step is based on decision which considers all previous data:

where . For simplicity, many works assume , that is, function evaluations are deterministic. However, we can easily extend the description to include stochastic functions (e.g.: homoscedastic noise ).

The search method is the sequence of decisions , which leads to the final decision . In most applications, the objective is to optimize the response of the final decisions. Then, the criteria relies on the optimality error or optimality gap, which can be expressed as:

In other applications, the objective may require to converge to in the input space. Then, we can use for example the Euclidean distance error:

The previous equations can also be interpreted as variants of the loss function for the decision at each step. Thus, the optimal decision is defined as the function that minimizes the loss function:

This requires full knowledge of function , which is unavailable. Instead, let assume that the target function belongs to a family of functions , e.g.: continuous functions in . Let also assume that the function can be represented as sample from a probability distribution over functions . Then, the best response case analysis for the search process is defined as the decision that optimizes the expectation of the loss function:

where is a prior distribution over functions.

However, we can improve the equation considering that, at decision we have already observed the actual response of the function at points, . Thus, the prior information of the function can be updated with the observations and the Bayes rule:

In fact, we can actually rewrite the equation to represent the updates sequentially:

Thus, the previous equation can be rewritten as:

This equation is the root of Bayesian optimization, where the Bayesian part comes from the fact that we are computing the expectation with respect to the posterior distribution, also called belief, over functions. Therefore, Bayesian optimization is a memory-based optimization algorithm.

As commented before, most of the theory of Bayesian optimization is related to deterministic functions, we consider also stochastic functions, that is, we assume there might be a random error in the function output. In fact, evaluations can produce different outputs if repeated. In that case, the target function is the expected output. Furthermore, in a recent paper by[5] it has been shown that, even for deterministic functions, it is better to assume certain error in the observation. The main reason being that, in practice, there might be some mismodelling errors which can lead to instability of the recursion if neglected.

# Bayesian optimization general model

In order to simplify the description, we are going to use a special case of Bayesian optimization model defined previously which corresponds to the most common application. In subsequent Sections we will introduce some generalizations for different applications.

Without loss of generality, consider the problem of finding the minimum of an unknown real valued function , where is a compact space, . Let be a prior distribution over functions represented as a stochastic process, for example, a Gaussian process , with inputs and an associate kernel or covariance function . Let also assume that the target function is a sample of the stochastic process .

In order to find the minimum, the algorithm has a maximum budget of evaluations of the target function . The purpose of the algorithm is to find optimal decisions that provide a better performance at the end.

One advantage of using Gaussian processes as a prior distributions over functions is that new observations of the target function can be easily used to update the distribution over functions. Furthermore, the posterior distribution is also a Gaussian process . Therefore, the posterior can be used as an informative prior for the next iteration in a recursive algorithm.

In a more general setting, many authors have suggested to modify the standard zero-mean Gaussian process for different variations that include semi-parametric models[9][7][11][19], use of hyperpriors on the parameters [15][2][8], Student t processes[6][22][28], etc.

We use a generalized linear model of the form:

where

The term means a non-parametric process, which can make reference to a Gaussian process or a Student's t process . In both cases, is the observation noise variance, sometimes called nugget, and it is problem specific. Many authors decide to fix this value when the function is deterministic, for example, a computer simulation. However, as cleverly pointed out in[5], there might be more reasons to include this term appart from being the observation noise, for example, to consider model inaccuracies.

This model has been presented in different ways depending on the field where it was used:

• As a generalized linear model with heteroscedastic perturbation .
• As a nonparametric process of the form .
• As a semiparametric model

# Models and functions

This library was originally developed for as part of a robotics research project[15][14], where a Gaussian process with hyperpriors on the mean and signal covariance parameters. Then, the metamodel was constructed using the Maximum a Posteriory (MAP) of the parameters. By that time, it only supported one kernel function, one mean function and one criterion.

However, the library now has grown to support many more surrogate models, with different distributions (Gaussian processes, Student's-t processes, etc.), with many kernels and mean functions. It also provides different criteria (even some combined criteria) so the library can be used to any problem involving some bounded optimization, stochastic bandits, active learning for regression, etc.

## Surrogate models

As seen in Section modopt this library implements only one general regression model. However, we can assign a set of priors on the parameters of the model , (the kernel hyperparameter will be discussed in Section learnker). Thus, the options are:

• "sGaussianProcess": a standard Gaussian process where the hyperparameters are known.
• "sGaussianProcessML": a standard Gaussian process where the hyperparameters are estimated directly from data using maximum likelihood estimates.
• "sGaussianProcessNormal": a Gaussian process with a Normal prior on the mean function parameters and known .
• "sStudentTProcessJef": in this case we use the Jeffreys prior for and . This is a kind of uninformative prior which is invariant to reparametrizations. Once we set a prior on the posterior becomes a Student's t Process.
• "sStudentTProcessNIG": in this case we standard conjugate priors, that is, a Normal prior on and a Inverse Gamma on . Therefore, the posterior is again a Student's t process.

Gaussian processes are a very general model that can achieve good performance with a reasonable computational cost. However, Student's t processes, thanks to the hierarchical structure of priors, provide an even more general setup for a minor extra cost. Furthermore, the Student's t distribution is robust to outliers and heavy tails in the data.

## Kernel (covariance) models

One of the critical components of Gaussian and Student's t processes is the definition of the kernel function, which defines the correlation between points in the input space. As a correlation function, the kernel must satisfy a set of properties (e.g.: being positive definite). All the kernel models available and its combinations satisfy the kernel restrictions.

The functions with "ISO" in their name are isotropic function, that is, they share a single set of parameters for all the dimensions of the input space.

The functions with "ARD" in their name use Automatic Relevance Determination, that is, they use independent parameters for every dimension of the problem. Therefore, they can be use to find the relevance of the each feature in the input space. In the limit, this can be used for feature selection.

### Atomic kernels

• "kConst": a simple constant function.
• "kLinear", "kLinearARD": a linear function.
• "kMaternISO1", "kMaternISO3","kMaternISO5","kMaternARD1","kMaternARD3","kMaternARD5": Matern kernel functions. The number divided by 2 represents the order of the function. See[20] for a description.
• "kPolyX": Polynomial kernel function. X is a number 1-6 which represents the exponent of the function.
• "kSEARD","kSEISO": Squared exponential kernel, also known as Gaussian kernel.
• "kRQISO": Rational quadratic kernel, also known as Student's t kernel.

### Binary kernels

This kernels allow to combine some of the previous kernels.

• "kSum": Sum of kernels.
• "kProd": Product of kernels.

Note that the binary kernels only admits two terms. However, we can combine them for more complex operations. For example if we write:

"kSum(kMaternISO3,kSum(kRQISO,kProd(kPoly4,kConst))"

it represents the expresion: Matern(3) + RationalQuadratic + C*Polynomial^4

In this case, the vector of parameters is splited from left to right: 1 for the Matern function, 2 for the RQ function, 2 for polynomial function and 1 for the constant. If the vector of parameters have more or less than 6 elements, the system complains.

## Parametric (mean) functions

Although the nonparametric process is able to model a large amount of funtions, we can model the expected value of the nonparametric process as a parametric function. This parametric model will help to capture large offsets and global trends.

The usage is analogous to the kernel functions.

• "mZero","mOne","mConst": constant functions. For simplicity and because they are largely used, we provide special cases f(x) = 0 and f(x) = 1.
• "mLinear": linear function.
• "mSum": binary function which can be used to combine other functions.

## Selection criteria

As discussed in Introduction to Bayesian Optimization, one of the critical aspects for Bayesian optimization is the decision (loss) function. Unfortunately, the functions described there are unavailable, because they assume knowledge of the optimal value . However, we can define proxies for those functions.

Some criteria, such as the expected improvement and the lower confidence bound admits an annealed version "cXXXa". In that version, the parameter that is used to trade off exploration and exploitation changes over time to priorize exploration at the begining and exploitation at the end.

Many criteria depends on the prediction function, which can be a Gaussian or a Student's t distribution, depending on the surrogate model. However, the library includes all the criteria for both distributions, and the system automatically selected the correct one.

### Atomic criteria

• "cEI","cBEI","cEIa": The most extended and reliable algorithm is the Expected Improvement algorithm[16]. In this case we provide the general version from[24] which includes an exponent to trade off exploration and exploitation "cEI". For an annealed version of the exploration/exploitation trade off, use "cEIa". Whe also includes a variation from[17] which add a bias or threshold to the improvement "cBEI".
• "cLCB", "cLCBa": Another popular algorithm is the Lower Confidence Bound (LCB), or UCB in case of maximization. Introduced by[cox1992statistical] as Sequential Design for Optimization (SDO). Analogously, "cLCBa" represents an annealed version of the exploration/exploitation trade off.
• "cMI": A generalized version of the LCB criterion which relies on the mutual information. See[4]
• "cPOI": Probability of improvement, by[12]
• "cExpReturn","cThompsonSampling","cOptimisticSampling": This criteria are related with the predicted return of the function. The first one is literally the expected return of the function (mean value). The second one is based on the Thompson sampling (drawing a random sample from the predicted distribution). Finally, the optimistic sampling takes the minimum of the other two (mean vs random).
• "cAopt": This is based on the A-optimality criteria. It is the predicted variance at the query point. Thus, this criteria is intended for exploration of the input space, not for optimization.
• "cDistance": This criteria adds a cost to a query point based on the distance with respect to the previous evaluation. Combined with other criteria functions, it might provide a more realistic setup for certain applications[13]

### Combined criteria

• "cSum","cProd": Sum and product of different criteria functions.
• "cHedge", "cHedgeRandom": Bandit based selection of the best criteria based on the GP-Hedge algorithm[8]. It automatically learns based on the behaviour of the criteria during the optimization process. The original version "cHedge" uses the maximum expected return as a reward for each criteria. We add a variant "cHedgeRandom" where the reward is defined in terms of Thompson sampling.

In this case, the combined criteria admits more that two functions. For example:

"cHedge(cSum(cEI,cDistance),cLCB,cPOI,cOptimisticSampling)"

## Methods for learning the kernel parameters

The posterior distribution of the model, which is necessary to compute the criterion function, cannot be computed in closed form if the kernel hyperparameters are unknown. Thus, we need a find to approximate this posterior distribution conditional on the kernel hyperparameters.

First, we need to consider if we are going to use a full Bayesian approach or an empirical Bayesian approach. The first one, computes the full posterior distribution by propagation of the uncertainty of each element and hyperparameter to the posterior. In this case, it can be done by discretization of the hyperparameter space or by using MCMC (not yet implemented). In theory, it is more precise but the computation burden can be orders of magnitude higher. The empirical approach on the other hand computes a point estimate of the hyperparameters based on some score function and use it as a "true" value. Although the uncertainty estimation in this case might not be as accurate as the full Bayesian, the computation overhead is minimal.

For the score function, we need to find the likelihood function of the observed data for the parameters. Depending on the model, this function will be a multivariate Gaussian distribution or multivariate t distribution. In general, we present the likelihood as a log-likelihood function up to a constant factor, that is, we remove the terms independent of from the log likelihood. In practice, whether we use a point estimate (maximum score) or full Bayes MCMC/discrete posterior, the constant factor is not needed.

We are going to consider the following score functions to learn the kernel hyperparameters:

• Leave one out cross-validation (SC_LOOCV): In this case, we try to maximize the average predicted log probability by the leave one out (LOO) cross-validation strategy. This is sometimes called a pseudo-likelihood.
• Maximum Total Likelihood (SC_MTL) For any of the models presented, one approach to learn the hyperparameters is to maximize the likelihood of all the parameters , and . Then, the likelihood function is a multivariate Gaussian distribution. We can obtain a better estimate if we adjust the number of degrees of freedom, this is called restricted maximum likelihood. The library automatically selects the restricted version, if it is suitable.
• Posterior maximum likelihood (SC_ML): In this case, the likelihood function is modified to consider the posterior estimate of based on the different cases defined in Section surrmods. In this case, the function will be a multivariate Gaussian or t distribution, depending on the kind of prior used for .
• Maximum a posteriori (SC_MAP): We can modify the previous algorithms by adding a prior distribution . By default, we add a joint normal prior on all the kernel hyperparameters. However, if the variance of the prior hp_std is invalid (<=0), then we assume a flat prior on that hyperparameter. Since we assume that the hyperparameters are independent, we can apply priors selectively only to a small set.

## Initial design methods

In order to build a suitable surrogate function, we a need a preliminar set of samples. In Bayesian optimization this is typically performed using alternative experimental design criteria. In this first step, usually the main criteria is space filling. Thus, we have implemented the subsequent designs:

• Uniform sampling: The search space is sampled uniformly.

Note: Since we do not assume any struture in the set of discrete points during discrete optimization, only uniform sampling of the discrete set is available in that case.


text:
# The Non-metricity Formulation of General Relativity

The Non-metricity Formulation of General Relativity After recalling the differential geometry of non-metric connections in the formalism of differential forms, we introduce the idea of a non-metricity (NM) connection, whose connection 1-forms coincides with the non-metricity 1-forms for a class of cobase fields. Then we formulate a theory of gravitation [(equivalent to General Relativity (GR)] which admits a geometrical interpretation in a flat torsionless space where the gravitational field is completely manifest in the non-metricity of a NM connection. We define and then apply the non-metricity gauge to a gravitational Lagrangian density discovered by Wallner (Acta Phys Austr 54:165–189, 1981) (proved in Appendix A to be equivalent to Einstein–Hilbert). The Einstein equations coupled to the matter currents $$\left( \mathcal {J}_{\alpha }\right)$$ J α thus becomes $$\delta dg_{\alpha }=\mathcal {T}_{\alpha }+\mathcal {J}_{\alpha }$$ δ d g α = T α + J α , where $$\left( \mathcal {T}_{\alpha }\right)$$ T α is identified as the gravitational energy-momentum currents, to which we shall find a relatively simple and physically appealing form. It is also shown that in the gravitational analogue of the Lorenz gauge, our field equations can be written as a system of Proca equations, which may be of interest in the study of propagation of gravitational-electromagnetic waves. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Advances in Applied Clifford Algebras Springer Journals

# The Non-metricity Formulation of General Relativity

, Volume 27 (3) – Jan 3, 2017
32 pages

/lp/springer_journal/the-non-metricity-formulation-of-general-relativity-b0uY6hcFIp
Publisher
Springer International Publishing
Subject
Physics; Mathematical Methods in Physics; Theoretical, Mathematical and Computational Physics; Applications of Mathematics; Physics, general
ISSN
0188-7009
eISSN
1661-4909
D.O.I.
10.1007/s00006-016-0749-8
Publisher site
See Article on Publisher Site

### Abstract

After recalling the differential geometry of non-metric connections in the formalism of differential forms, we introduce the idea of a non-metricity (NM) connection, whose connection 1-forms coincides with the non-metricity 1-forms for a class of cobase fields. Then we formulate a theory of gravitation [(equivalent to General Relativity (GR)] which admits a geometrical interpretation in a flat torsionless space where the gravitational field is completely manifest in the non-metricity of a NM connection. We define and then apply the non-metricity gauge to a gravitational Lagrangian density discovered by Wallner (Acta Phys Austr 54:165–189, 1981) (proved in Appendix A to be equivalent to Einstein–Hilbert). The Einstein equations coupled to the matter currents $$\left( \mathcal {J}_{\alpha }\right)$$ J α thus becomes $$\delta dg_{\alpha }=\mathcal {T}_{\alpha }+\mathcal {J}_{\alpha }$$ δ d g α = T α + J α , where $$\left( \mathcal {T}_{\alpha }\right)$$ T α is identified as the gravitational energy-momentum currents, to which we shall find a relatively simple and physically appealing form. It is also shown that in the gravitational analogue of the Lorenz gauge, our field equations can be written as a system of Proca equations, which may be of interest in the study of propagation of gravitational-electromagnetic waves.

### Journal

Advances in Applied Clifford AlgebrasSpringer Journals

Published: Jan 3, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations

Abstract access only

18 million full-text articles

Print

20 pages / month

PDF Discount

20% off

format: plain


text:
# Install issue: Mathematica 9 on Ubuntu 13.10

I am trying to install Mathematica 9 on a Ubuntu 13.10 system. I had this working on my previous 13.04 version but now it breaks on when trying to unpack the sh archive:

\$ sh Mathematica_9.0.0_LINUX.sh
Mathematica Secured 9.0.0 for LINUX Installer Archive

Verifying archive integrity.
Extracting installer. ...........................................................................
CRITICAL FAILURE: Fundamental Error

sudo sh ... produces same error, the md5 hash checks out ok.

Anybody encountered this problem, and more to the point, figured out what is wrong? I submitted a ticket to Wolfram support but no luck there so far.

Thanks for any ideas to try.

-
Make sure there is no space in your directory name. –  alephalpha Feb 19 at 5:19
Yes I tried making it executable (should't make a difference and it didn't). Yes 'sh' is the right shell for interpreting the archive. Sticking point now is that my activation code doesn't work since the fingerprint of the machine is different after OS upgrade... –  SEngstrom Feb 19 at 5:20
@alephalpha as you can see there are no spaces in the directory I was working from. –  SEngstrom Feb 19 at 5:21
try renaming the .sh file to mma9.sh (something without dot and dash) and then run it. –  Hubble07 Feb 19 at 9:38

I worked around this by directing the shell archive extraction to keep the files:

sh Mathematica_9.0.0_LINUX.sh --keep

Then I cd'd into the Installation directory and ran the installer from there. Note that the temporary directory will named after the process number you ran the extraction in so will differ from what I had:

cd .4725/Unix/Installer
sudo MathInstaller

... success. Now to convince Wolfram to give me a new activation code for the upgraded machine.

-
Now upgrading to Mathematica 10.0.0 - Archive extraction still breaks, this time under Ubuntu 14.04 LTS. Fix above still seems to work. –  SEngstrom Jul 11 at 16:19


text:
# How to type superscript correctly

Hi there,

Thanks for helping me out with my language problem@Erry_jonson

I thought it is exceptionally simple yet I was unable to make it.

I was attempting to type cm-2(-2 as superscript) on my item depiction. It works on Microsoft Word yet on my site, it became “cm-2” eventually.

Then, at that point, I attempted to type thusly: “cm^-2”. Subsequent to distributing it became “cm￪2” some way or another the “^” and "- " have been joined together to a bolt…

Eventually, I composed “cm^ - 2”, definitely I composed a “space” between

there are my inquiries:

how is it possible that I would type “- 2” as a superscript? (not 2)

how should web-flow not naturally turn “^-” into “￪”?

Thank you.
Dij

hi @Dij_L you are misunderstand what is and how content of websites works. I will recommend to invest some time to learn what is HTML and CSS and how these technologies differ from standard software for writing text (text editors).

Anyway back to you question if you browse internet with keywords  html superscript or how to superscript website you will find many great sources that will reveal secret weapon of  .

if you will write more complex math formula on website there are libraries that can may help you.

http://mathquill.com

but feel free to fine one that will fulfil your expectations

because HTML is not text editor software

1 Like

Hey @Dij_L , on the note from @Stan ,  support is not available on Webflow as of now. A workaround I tried may work for you:

1. Write out cm-2 normally
2. Wrap -2 within a span and give it a class superscript
3. Style this class by reducing font-size to 10 or your desired size, change display to to inline-block.
4. Switching display to inline-block will allow css transform, so you can create a transform move on y-axis.

Let me know if this works!

2 Likes

hi @imtiazraqib Im sorry that  doesn’t work for you so you can also use ::after approach. But before you try these complicated ways to create superscript just use .

let me know if this works for you !!!

1 Like

Hi @Stan and @imtiazraqib, thanks for all the tips.
Actually this subscript belongs to one paragraph of a dynamic product description. So at beginning I was expecting to type superscript without HTML Embed. For my problem, I could type 15 N/cm² instead of 15 N cm-2. But I am curious why I could just type cm² but not cm-²(after typed “-”, “cm²” became “cm-2”).

I tried your solutions, all work on a normal page but I am still trying to make it works with dynamic content.

thank you again, I will keep update
Cheers

hi dude
that worked for me =)
I do not know why it did not work for guys
But it worked for me :)))
thanks

1 Like

hi @Dij_L  is working in CMS RTE as well. Here is just screenshot to prove it. I don’t know what other way to prove that  is working and if you don’t know how to use it just use @imtiazraqib approach or try to use another approach with use of ::after pseudo class.

2 Likes

HI @Dij_L if you don’t have further question related to this topic feel free to close your request as solved

This topic was automatically closed 24 hours after the last reply. New replies are no longer allowed.


text:
# Math Help - Indepedant events.

1. ## Independant events.

I hope this isn't too long winded, I have most the answers, just need a check and some help on the last part.

Laser surgery to fix short-sightedness is becoming more popular. However, for some people, a second procedure is necessary. The following table lists the joint probability of needing a second procedure and whether the patient has a corrective lens with a factor (dioptre) of minus 8 or less.

This is supposed to be a table!

__________________ Vision corrective factor > -8, Vision corrective factor <= -8

1st procedure successful $~~~~~~$___________0.66 $~~~~~~$___________________0.15

2nd procedure required $~~~~~~$____________0.05 $~~~~~~$___________________0.14

a) Find the probability that a second procedure is required.

This is simply $0.05+0.14 = 0.19$

b) Determine the probability that someone whose corrective lens factor is minus 8 or less does not require a second procedure.

This is $0.66+0.15+0.14 = 0.95$

I think the answer is yes, but I don't know how to explain

2. Originally Posted by Bushy
I hope this isn't too long winded, I have most the answers, just need a check and some help on the last part.

Laser surgery to fix short-sightedness is becoming more popular. However, for some people, a second procedure is necessary. The following table lists the joint probability of needing a second procedure and whether the patient has a corrective lens with a factor (dioptre) of minus 8 or less.

This is supposed to be a table!

__________________ Vision corrective factor > -8, Vision corrective factor <= -8

1st procedure successful $~~~~~~$___________0.66 $~~~~~~$___________________0.15

2nd procedure required $~~~~~~$____________0.05 $~~~~~~$___________________0.14

a) Find the probability that a second procedure is required.

This is simply $0.05+0.14 = 0.19$

b) Determine the probability that someone whose corrective lens factor is minus 8 or less does not require a second procedure.

This is $0.66+0.15+0.14 = 0.95$


text:
Feeds:
Posts

## 08. An arithmetic partition

08. PEN O 35 ( Romania TST 1998 )

Let $n \ge 3$ be a prime number and $a_{1} < a_{2} < \cdots < a_{n}$ be integers.
Prove that $a_{1}, \cdots,a_{n}$ is an arithmetic progression if and only if there exists a partition of $\{0, 1, 2, \cdots \}$ into sets $A_{1},A_{2},\cdots,A_{n}$ such that

$a_{1} + A_{1} = a_{2} + A_{2} = \cdots = a_{n} + A_{n}$

where $x + A$ denotes the set $\{x + a \vert a \in A \}$.


text:
## What Is Stochastic Rounding?

In finite precision arithmetic the result of an elementary arithmetic operation does not generally lie in the underlying number system, $F$, so it must be mapped back into $F$ by the process called rounding. The most common choice is round to nearest, whereby the nearest number in $F$ to the given number is chosen; this is the default in IEEE standard floating-point arithmetic.

Round to nearest is deterministic: given the same number it always produces the same result. A form of rounding that randomly rounds to the next larger or next smaller number was proposed Barnes, Cooke-Yarborough, and Thomas (1951), Forysthe (1959), and Hull and Swenson (1966). Now called stochastic rounding, it comes in two forms. The first form rounds up or down with equal probability $1/2$.

We focus here on the second form of stochastic rounding, which is more interesting. To describe it we let $x\notin F$ be the given number and let $x_1$ and $x_2$ with $x_1 < x_2$ be the adjacent numbers in $F$ that are the candidates for the result of rounding. The following diagram shows the situation where $x$ lies just beyond the half-way point between $x_1$ and $x_2$:

We round up to $x_2$ with probability $(x-x_1)/(x_2-x_1)$ and down to $x_1$ with probability $(x_2-x)/(x_2-x_1)$; note that these probabilities sum to $1$. The probability of rounding to $x_1$ and $x_2$ is proportional to one minus the relative distance from $x$ to each number.

If stochastic rounding chooses to round to $x_1$ in the diagram then it has committed a larger error than round to nearest, which rounds to $x_2$. Indeed if we focus on floating-point arithmetic then the result of the rounding is

$\mathrm{f\kern.2ptl}(x) = x (1+\delta), \quad |\delta|\le 2u,$

where $u$ is the unit roundoff. For round to nearest the same result is true with the smaller bound $|\delta|\le u$. In addition, certain results that hold for round to nearest, such as $\mathrm{f\kern.2ptl}(x) = -\mathrm{f\kern.2ptl}(-x)$ and $\mathrm{f\kern.2ptl}(x/\sqrt{x^2+y^2}) \le 1$, can fail for stochastic rounding. What, then, is the benefit of stochastic rounding?

It is not hard to show that the expected value of the result of stochastically rounding $x$ is $x$ itself (this is obvious if $x = (x_1+x_2)/2$, for example), so the expected error is zero. Hence stochastic rounding maintains, in a statistical sense, some of the information that is discarded by a deterministic rounding scheme. This property leads to some important benefits, as we now explain.

In certain situations, round to nearest produces correlated rounding errors that cause systematic error growth. This can happen when we form the inner product of two long vectors $x$ and $y$ of nonnegative elements. If the elements all lie on $[0,1]$ then clearly the partial sum can grow monotonically as more and more terms are accumulated until at some point all the remaining terms “drop off the end” of the computed sum and do not change it—the sum stagnates. This phenomenon was observed and analyzed by Huskey and Hartree as long ago as 1949 in solving differential equations on the ENIAC. Stochastic rounding avoids stagnation by rounding up rather than down some of the time. The next figure gives an example. Here, $x$ and $y$ are $n$-vectors sampled from the uniform $[0,1]$ distribution, the inner product is computed in IEEE standard half precision floating-point arithmetic ($u \approx 4.9\times 10^{-4}$), and the backward errors are plotted for increasing $n$. From about $n = 1000$, the errors for stochastic rounding (SR, in orange) are smaller than those for round to nearest (RTN, in blue), the latter quickly reaching 1.

The figure is explained by recent results of Connolly, Higham, and Mary (2020). The key property is that the rounding errors generated by stochastic rounding are mean independent (a weaker version of independence). As a consequence, an error bound proportional to $\sqrt{n}u$ can be shown to hold with high probability. With round to nearest we have only the usual worst-case error bound, which is proportional to $nu$. In the figure above, the solid black line is the standard backward error bound $nu$ while the dashed black line is $\sqrt{n}u$. In this example the error with round to nearest is not bounded by a multiple of $\sqrt{n}u$, as the blue curve crosses the dashed black one. The underlying reason is that the rounding errors for round to nearest are correlated (they are all negative beyond a certain point in the sum), whereas those for stochastic rounding are uncorrelated.

Another notable property of stochastic rounding is that the expected value of the computed inner product is equal to the exact inner product. The same is true more generally for matrix–vector and matrix–matrix products and the solution of triangular systems.

Stochastic rounding is proving popular in neural network training and inference when low precision arithmetic is used, because it avoids stagnation.

Implementing stochastic rounding efficiently is not straightforward, as the definition involves both a random number and an accurate value of the result that we are trying to compute. Possibilities include modifying low level arithmetic routines (Hopkins et al., 2020) and exploiting the augmented operations in the latest version of the IEEE standard floating-point arithmetic (Fasi and Mikaitis, 2020).

Hardware is starting to become available that supports stochastic rounding, including the Intel Lohi neuromorphic chip, the Graphcore Intelligence Processing Unit (intended to accelerate machine learning), and the SpiNNaker2 chip.

Stochastic rounding can be done in MATLAB using the chop function written by me and Srikara Pranesh. This function is intended for experimentation rather than efficient computation.

## References

This is a minimal set of references, which contain further useful references within.

## What Is the Hilbert Matrix?

The Hilbert matrix $H_n = (h_{ij})$ is the $n\times n$ matrix with $h_{ij} = 1/(i+j-1)$. For example,

$H_4 = \left[\begin{array}{@{\mskip 5mu}c*{3}{@{\mskip 15mu} c}@{\mskip 5mu}} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\[6pt] \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5}\\[6pt] \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6}\\[6pt] \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7}\\[6pt] \end{array}\right].$

It is probably the most famous test matrix and its conditioning and other properties were extensively studied in the early days of digital computing, especially by John Todd.

The Hilbert matrix is symmetric and it is a Hankel matrix (constant along the anti-diagonals). Less obviously, it is symmetric positive definite (all its eigenvalues are positive) and totally positive (every submatrix has positive determinant). Its condition number grows rapidly with $n$; indeed for the 2-norm the asymptotic growth rate is $\kappa_2(H_n) \sim \mathrm{e}^{3.5n}$. An underlying reason for the ill conditioning is that the Hilbert matrix is obtained when least squares polynomial approximation is done using the monomial basis, and this is known to be an ill-conditioned polynomial basis.

The inverse of $H_n$ has integer entries, with $(i,j)$ entry

$(-1)^{i+j} (i+j-1) \displaystyle{ n+i-1 \choose n-j } { n+j-1 \choose n-i } { i+j-2 \choose i-1 }^2.$

For example,

$H_4^{-1} = \left[\begin{array}{rrrr} 16 & -120 & 240 & -140 \\ -120 & 1200 & -2700 & 1680 \\ 240 & -2700 & 6480 & -4200 \\ -140 & 1680 & -4200 & 2800 \\ \end{array}\right].$

The alternating sign pattern is a consequence of total positivity.

The determinant is also known explicitly:

$\det(H_n) = \displaystyle\frac{ (1!\, 2!\, \dots (n-1)!\, )^4 }{ 1!\, 2!\, \dots (2n-1)! } \sim 2^{-2n^2}. % \cite{todd54}$

The Hilbert matrix is infinitely divisible, which means that the matrix with $(i,j)$ element $h_{ij}^t$ is positive semidefinite for all nonnegative real numbers $t$.

Other interesting matrices can be formed from the Hilbert matrix. The Lotkin matrix, A = gallery('lotkin',n) in MATLAB, is the Hilbert matrix with the first row replaced by $1$s, making it nonsymmetric. Like the Hilbert matrix, it is ill conditioned and has an inverse with integer entries. For the $n\times n$ matrix with $(i,j)$ entry $1/(i+j)$, which is a submatrix of $H_{n+1}$, the largest singular value tends to $\pi$ as $n\to\infty$, with error decaying slowly as $O(1/\log n)$, as was shown by Taussky.

The Hilbert matrix is not a good test matrix, for several reasons. First, it is too ill conditioned, being numerically singular in IEEE double precision arithmetic for $n = 12$. Second, it is too special: the definiteness and total positivity, together with the fact that it is a special case of a Cauchy matrix (which has elements $1/(x_i + y_j)$), mean that certain quantities associated with it can be calculated more accurately than for a general positive definite matrix. The third reason why the Hilbert matrix is not a good test matrix is that most of its elements cannot be stored exactly in floating-point arithmetic, so the matrix actually stored is a perturbation of $H_n$—and since $H_n$ is ill conditioned this means (for example) that the inverse of the stored matrix can be very different from the exact inverse. A way around this might be to work the integer matrix $H_n^{-1}$, but its elements are exactly representable in double precision only for $n\le 12$.

MATLAB has functions hilb and invhilb for the Hilbert matrix and its inverse. How to efficiently form the Hilbert matrix in MATLAB is an interesting question. The hilb function essentially uses the code

j = 1:n;
H = 1./(j'+j-1)

This code is more concise and efficient than two nested loops would be. The second line may appear to be syntactically incorrect, because j'+j adds a column vector to a row vector. However, the MATLAB implicit expansion feature expands j' into an $n\times n$ matrix in which every column is j', expands j into an $n\times n$ matrix in which every row is j, then adds the two matrices. It is not hard to see that the Hilbert matrix is the result of these two lines of code.

## References

This is a minimal set of references, which contain further useful references within.

## Related Blog Posts

Posted in what-is | 3 Comments

## What Is a Fréchet Derivative?

Let $U$ and $V$ be Banach spaces (complete normed vector spaces). The Fréchet derivative of a function $f:U \to V$ at $X\in U$ is a linear mapping $L:U\to V$ such that

$f(X+E) - f(X) - L(X,E) = o(\|E\|)$

for all $E\in U$. The notation $L(X,E)$ should be read as “the Fréchet derivative of $f$ at $X$ in the direction $E$”. The Fréchet derivative may not exist, but if it does exist then it is unique. When $U = V = \mathbb{R}$, the Fréchet derivative is just the usual derivative of a scalar function: $L(x,e) = f'(x)e$.

As a simple example, consider $U = V = \mathbb{R}^{n\times n}$ and $f(X) = X^2$. From the expansion

$f(X+E) - f(X) = XE + EX + E^2$

we deduce that $L(X,E) = XE + EX$, the first order part of the expansion. If $X$ commutes with $E$ then $L_X(E) = 2XE = 2EX$.

More generally, it can be shown that if $f$ has the power series expansion $f(x) = \sum_{i=0}^\infty a_i x^i$ with radius of convergence $r$ then for $X,E\in\mathbb{R}^{n\times n}$ with $\|X\| < r$, the Fréchet derivative is

$L(X,E) = \displaystyle\sum_{i=1}^\infty a_i \displaystyle\sum_{j=1}^i X^{j-1} E X^{i-j}.$

An explicit formula for the Fréchet derivative of the matrix exponential, $f(A) = \mathrm{e}^A$, is

$L(A,E) = \displaystyle\int_0^1 \mathrm{e}^{A(1-s)} E \mathrm{e}^{As} \, ds.$

Like the scalar derivative, the Fréchet derivative satisfies sum and product rules: if $g$ and $h$ are Fréchet differentiable at $A$ then

\begin{alignedat}{2} f &= \alpha g + \beta h &&\;\Rightarrow\; L_f(A,E) = \alpha L_g(A,E) + \beta L_h(A,E),\\ f &= gh &&\;\Rightarrow\; L_f(A,E) = L_g(A,E) h(A) + g(A) L_h(A,E). \end{alignedat}

A key requirement of the definition of Fréchet derivative is that $L(X,E)$ must satisfy the defining equation for all $E$. This is what makes the Fréchet derivative different from the Gâteaux derivative (or directional derivative), which is the mapping $G:U \to V$ given by

$G(X,E) = \lim_{t\to0} \displaystyle\frac{f(X+t E)-f(X)}{t} = \frac{\mathrm{d}}{\mathrm{dt}}\Bigm|_{t=0} f(X + tE).$

Here, the limit only needs to exist in the particular direction $E$. If the Fréchet derivative exists at $X$ then it is equal to the Gâteaux derivative, but the converse is not true.

A natural definition of condition number of $f$ is

$\mathrm{cond}(f,X) = \lim_{\epsilon\to0} \displaystyle\sup_{\|\Delta X\| \le \epsilon \|X\|} \displaystyle\frac{\|f(X+\Delta X) - f(X)\|}{\epsilon\|f(X)\|},$

and it can be shown that $\mathrm{cond}$ is given in terms of the Fréchet derivative by

$\mathrm{cond}(f,X) = \displaystyle\frac{\|L(X)\| \|X\|}{\|f(X)\|},$

where

$\|L(X)\| = \sup_{Z\ne0}\displaystyle\frac{\|L(X,Z)\|}{\|Z\|}.$

For matrix functions, the Fréchet derivative has a number of interesting properties, one of which is that the eigenvalues of $L(X)$ are the divided differences

$f[\lambda_i,\lambda_j] = \begin{cases} \dfrac{ f(\lambda_i)-f(\lambda_j) }{\lambda_i - \lambda_j}, & \lambda_i\ne\lambda_j, \\ f'(\lambda_i), & \lambda_i=\lambda_j, \end{cases}$

for $1\le i,j \le n$, where the $\lambda_i$ are the eigenvalues of $X$. We can check this formula in the case $F(X) = X^2$. Let $(\lambda,u)$ be an eigenpair of $X$ and $(\mu,v)$ an eigenpair of $X^T$, so that $Xu = \lambda u$ and $X^Tv = \mu v$, and let $E = uv^T$. Then

$L(X,E) = XE + EX = Xuv^T + uv^TX = (\lambda + \mu) uv^T.$

So $uv^T$ is an eigenvector of $L(X)$ with eigenvalue $\lambda + \mu$. But $f[\lambda,\mu] = (\lambda^2-\mu^2)/(\lambda - \mu) = \lambda + \mu$ (whether or not $\lambda$ and $\mu$ are distinct).

## References

This is a minimal set of references, which contain further useful references within.

## Six Useful LaTeX Packages

Here are six $\LaTeX$ packages that, while they are perhaps not among the best known, I have found to be very useful. All of them are probably available in your $\LaTeX$ distribution, but in case they are not I give links to the packages on CTAN (the Comprehensive TeX Archive Network). A PDF file containing a demonstration of the examples below is available here.

## Mathtools

Mathtools extends the amsmath package with various additions, particularly for fine-tuning the typesetting of equations, and it also fixes a few problems with amsmath. To give just one example, the standard amsmath bmatrix environment centres its columns:

A = \begin{bmatrix}
12 & 3  & 0 \\
0  & 10 &-1
\end{bmatrix}

For matrices with negative elements, such as this one, the columns can look better right justified. The bmatrix* environment provided by mathtools allows a column alignment parameter to be given, so we can right justify the columns as follows:

A = \begin{bmatrix*}[r]
12 & 3  & 0 \\
0  & 10 &-1
\end{bmatrix*}

Unfortunately, only one alignment parameter can be given, so this command does not allow different columns to have different alignments (to achieve this, we can use the array environment instead).

The mathtools package is regularly updated.

## Url

The url package provides the \url command, which acts like \verb except that it allows line breaks to be made when the argument is typeset. I use url whenever I need to typeset a URL or a path, in particular in my BibTeX bib files. I call the package with the hyphens option, so that line breaks will be allowed at hyphens. Example:

\usepackage[hyphens]{url}
...
\footnote{\url{https://icl.bitbucket.io/hpl-ai/}}
\url{http://blogs.mathworks.com/cleve/2013/01/07/george-forsythe/}

An advantage over the \verb command is that \url (or, if necessary, \urldef) can be used inside another command, such as \footnote, as in this example.

## Upref

By loading the upref package you ensure that every \ref or \eqref reference is set in an upright font. This is normal typesetting practice, but it is not enforced by all $\LaTeX$ classes. It is enforced by the SIAM class, for example.

## Siunitx

The siunitx package typesets SI units according to the official rules, with just the right spacing between the number and the unit. For example, you can type

\si{kg.m/s^2}

$\mathrm{kg}\,\mathrm{m}/\mathrm{s}^2$

Even if you do not need to type SI units, you might find useful the package’s ability to round and format numbers.

## Upquote

The upquote package forces left and right quotes to be typed as straight (vertical) quotes in \verb and verbatim environments. This is usually what is required in program listings. We used upquote in MATLAB Guide.

## Fancyvrb

The fancyvrb package provides a more powerful replacement for the verbatim environment, allowing the font to be specified, the environment to be boxed, an external file to be inclued in a verbatim environment, and much more. Here is example that formats the contents of the environment in blue and boldface:

\begin{Verbatim}[fontsize=\small,formatcom=\color{blue}]
>> syms x, int(1/(1-x^4))
ans =
atan(x)/2 + atanh(x)/2
\end{Verbatim}

The second example includes an external file, numbers the lines, and puts horizontal lines before and after the code:

\renewcommand{\theFancyVerbLine}%
{\textcolor{red}{\small\arabic{FancyVerbLine}}}
\VerbatimInput[frame=lines,numbers=left,fontsize=\small]%
{d:/matlab/gen/gf.m}

## Other packages

I have previously written about the enumitem and booktabs packages, which I also strongly recommend.

## What Is the Adjugate of a Matrix?

The adjugate of an $n\times n$ matrix $A$ is defined by

$\mathrm{adj}(A) = \bigl( (-1)^{i+j} \det(A_{ji}) \bigr),$

where $A_{pq}$ denotes the submatrix of $A$ obtained by deleting row $p$ and column $q$. It is the transposed matrix of cofactors. The adjugate is sometimes called the (classical) adjoint and is sometimes written as $A^\mathrm{A}$.

For nonsingular $A$,

$\mathrm{adj}(A) = \det(A) A^{-1},$

and this equation is often used to give a formula for $A^{-1}$ in terms of determinants.

Since the adjugate is a scalar multiple of the inverse, we would expect it to share some properties of the inverse. Indeed we have

\begin{aligned} \mathrm{adj}(AB) &= \mathrm{adj}(B) \mathrm{adj}(A),\\ \mathrm{adj}(A^*) &= \mathrm{adj}(A)^*. \end{aligned}

These properties can be proved by first assuming the matrices are nonsingular—in which case they follow from properties of the determinant and the inverse—and then using continuity of the entries of $\mathrm{adj}(A)$ and the fact that every matrix is the limit of a sequence of nonsingular matrices. Another property that can be proved in a similar way is

$\mathrm{adj}\bigl(\mathrm{adj}(A)\bigr) = (\det A)^{n-2} A.$

If $A$ has rank $n-2$ or less then $\mathrm{adj}(A)$ is just the zero matrix. Indeed in this case every $(n-1)\times(n-1)$ submatrix of $A$ is singular, so $\det(A_{ji}) \equiv 0$ in the definition of $\mathrm{adj}(A)$. In particular, $\mathrm{adj}(0) = 0$ and $\mathrm{adj}(xy^*) = 0$ for any rank-1 matrix $xy^*$

Less obviously, if $\mathrm{rank}(A) = n-1$ then $\mathrm{adj}(A)$ has rank 1. This can be seen from the formula below that expresses the adjugate in terms of the singular value decomposition (SVD).

If $x,y\in\mathbb{C}^n$, then for nonsingular $A$,

\begin{aligned} \det(A + xy^*) &= \det\bigl( A(I + A^{-1}x y^*) \bigr)\\ &= \det(A) \det( I + A^{-1}x y^*)\\ &= \det(A) (1+ y^*A^{-1}x)\\ &= \det(A) + y^* (\det(A)A^{-1})x\\ &= \det(A) + y^* \mathrm{adj}(A)x. \end{aligned}

Again it follows by continuity that for any $A$,

$\det(A + xy^*) = \det(A) + y^* \mathrm{adj}(A)x.$

This expression is useful when we need to deal with a rank-1 perturbation to a singular matrix. A related identity is

$\det\left( \begin{bmatrix} A & x \\ y^* & \alpha \\ \end{bmatrix}\right) = \alpha \det(A) - y^* \mathrm{adj}(A) x.$

Let $A = U\Sigma V^*$ be an SVD, where $U$ and $V$ are unitary and $\Sigma = \mathrm{diag}(\sigma_i)$, with $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0$. Then

$\mathrm{adj}(A) = \mathrm{adj}(V^*) \mathrm{adj}(\Sigma) \mathrm{adj}(U) = \mathrm{adj}(V)^* \mathrm{adj}(\Sigma) \mathrm{adj}(U).$

It is easy to see that $D = \mathrm{adj}(\Sigma)$ is diagonal, with

$d_{kk} = \displaystyle\prod_{i=1\atop i\ne k }^n \sigma_i.$

Since $U$ and $V$ are unitary and hence nonsingular,

\begin{aligned} \mathrm{adj}(A) &= \mathrm{adj}(V)^* \mathrm{adj}(\Sigma) \mathrm{adj}(U)\\ &= \overline{\det(V)} V^{-*} \mathrm{adj}(\Sigma) \det(U) U^{-1}\\ &= \overline{\det(V)} V \mathrm{adj}(\Sigma) \det(U) U^*\\ &= \overline{\det(V)} \det(U) \cdot V \mathrm{adj}(\Sigma) U^*. \end{aligned}

If $\mathrm{rank}(A) = n-1$ then $\sigma_n = 0$ and so

$\mathrm{adj}(A) = \rho\, \sigma_1 \sigma_2 \cdots \sigma_{n-1}v_n^{} u_n^*,$

where $\rho = \overline{\det(V)} \det(U)$ has modulus $1$ and $v_n$ and $u_n$ are the last columns of $V$ and $U$, respectively.

The definition of $\mathrm{adj}(A)$ does not provide a good means for computing it, because the determinant computations are prohibitively expensive. The following MATLAB function (which is very similar to the function adjoint in the Symbolic Math Toolbox) uses the SVD formula.

function X = adj(A)
%   X = ADJ(A) computes the adjugate of the square matrix A via
%   the singular value decomposition.

n = length(A);
[U,S,V] = svd(A);
D = zeros(n);
for i=1:n
d = diag(S);
d(i) = 1;
D(i,i) = prod(d);
end
X = conj(det(V))*det(U)*V*D*U';

Note that, in common with most SVD codes, the svd function does not return $\det(U)$ and $\det(V)$, so we must compute them. This function is numerically stable, as shown by Stewart (1998).

Finally, we note that Richter (1954) and Mirsky (1956) obtained the Frobenius norm bound

$\|\mathrm{adj}(A) \|_F \le \displaystyle\frac{ \|A\|_F^{n-1} }{ n^{(n-2)/2} }.$

## References

This is a minimal set of references, which contain further useful references within.

Posted in what-is | 2 Comments

## What Is a Matrix Function?

If you give an array as the argument to a mathematical function in a programming language or problem solving environment you are likely to receive the result of applying the function to each component of the array. For example, here MATLAB exponentiates the integers from 0 to 3:

>> A = [0 1; 2 3], X = exp(A)
A =
0     1
2     3
X =
1.0000e+00   2.7183e+00
7.3891e+00   2.0086e+01

When the array represents a matrix this componentwise evaluation is not useful, as it does not obey the rules of matrix algebra. To see what properties we might require, consider the matrix square root. This function is notable historically because Cayley considered it in the 1858 paper in which he introduced matrix algebra.

A square root of an $n\times n$ matrix $A$ is a matrix $X$ such that $X^2 = A$. Any square root $X$ satisfies

$AX = X^2 X = X X^2 = XA,$

so $X$ commutes with $A$. Also, $A^* = (X^2)^* = (X^*)^2$, so $X^*$ is a square root of $A^*$ (here, the superscript $*$ denotes the conjugate transpose).

Furthermore, for any nonsingular matrix $Z$ we have

$(Z^{-1}XZ)^2 = Z^{-1}X^2 Z = Z^{-1}A Z.$

If we choose $Z$ as a matrix that takes $A$ to its Jordan canonical form $J$ then we have $(Z^{-1}XZ)^2 = J$, so that $Z^{-1}XZ$ is a square root of $J$, or in other words $X$ can be expressed as $X = Z \sqrt{J} Z^{-1}$.

Generalizing from these properties of the matrix square root it is reasonable to ask that a function $f(A)$ of a square matrix $A$ satisfies

• $f(A)$ commutes with $A$,
• $f(A^*) = f(A)^*$,
• $f(A) = Z f(J)Z^{-1}$, where $A$ has the Jordan canonical form $A = ZJZ^{-1}$.

(These are of course not the only requirements; indeed $f(A) \equiv I$ satisfies all three conditions.)

Many definitions of $f(A)$ can be given that satisfy these and other natural conditions. We will describe just three definitions (a notable omission is a definition in terms of polynomial interpolation).

## Power Series

For a function $f$ that has a power series expansion we can define $f(A)$ by substituting $A$ for the variable. It can be shown that the matrix series will be convergent for matrices whose eigenvalues lie within the radius of convergence of the scalar power series. For example, where $\rho$ denotes the spectral radius. This definition is natural for functions that have a power series expansion, but it is rather limited in its applicability.

\begin{aligned} \cos(A) &= I - \displaystyle\frac{A^2}{2!} + \frac{A^4}{4!} - \frac{A^6}{6!} + \cdots,\\ \log(I+A) &= A - \displaystyle\frac{A^2}{2} + \frac{A^3}{3} - \frac{A^4}{4} + \cdots, \quad \rho(A)<1, \end{aligned}

## Jordan Canonical Form Definition

If $A$ has the Jordan canonical form

$Z^{-1}AZ = J = \mathrm{diag}(J_1, J_2, \dots, J_p),$

where

$J_k = J_k(\lambda_k) = \begin{bmatrix} \lambda_k & 1 & & \\ & \lambda_k & \ddots & \\ & & \ddots & 1 \\ & & & \lambda_k \end{bmatrix} = \lambda_k I + N \in \mathbb{C}^{m_k\times m_k},$

then

$f(A) := Z f(J) Z^{-1} = Z \mathrm{diag}(f(J_k)) Z^{-1},$

where

$f(J_k) := \begin{bmatrix} f(\lambda_k) & f'(\lambda_k) & \dots & \displaystyle\frac{f^{(m_k-1)}(\lambda_k)}{(m_k-1)!} \\ & f(\lambda_k) & \ddots & \vdots \\ & & \ddots & f'(\lambda_k) \\ & & & f(\lambda_k) \end{bmatrix}.$

Notice that $f(J_k)$ has the derivatives of $f$ along its diagonals in the upper triangle. Write $J_k = \lambda_k I + N$, where $N$ is zero apart from a superdiagonal of 1s. The formula for $f(J_k)$ is just the Taylor series expansion

$f(J_k) = f(\lambda_k I + N) = f(\lambda_k)I + f'(\lambda_k)N + \displaystyle\frac{f''(\lambda_k)}{2!}N^2 + \cdots + \displaystyle\frac{f^{(m_k-1)}(\lambda_k)}{(m_k-1)!}N^{m_k-1}.$

The series is finite because $N^{m_k} = 0$: as $N$ is powered up the superdiagonal of 1s moves towards the right-hand corner until it eventually disappears.

An immediate consequence of this formula is that $f(A)$ is defined only if the necessary derivatives exist: for each eigenvalue $\lambda_k$ we need the existence of the derivatives at $\lambda_k$ of order up to one less than the size of the largest Jordan block in which $\lambda_k$ appears.

When $A$ is a diagonalizable matrix the definition simplifies considerably: if $A = ZDZ^{-1}$ with $D = \mathrm{diag}(\lambda_i)$ then $f(A) = Zf(D)Z^{-1} = Z \mathrm{diag}(\lambda_i) Z^{-1}$.

## Cauchy Integral Formula

For a function $f$ analytic on and inside a closed contour $\Gamma$ that encloses the spectrum of $A$ the matrix $f(A)$ can be defined by the Cauchy integral formula

$f(A) := \displaystyle\frac{1}{2\pi \mathrm{i}} \int_\Gamma f(z) (zI-A)^{-1} \,\mathrm{d}z.$

This formula is mathematically elegant and can be used to provide short proofs of certain theoretical results.

## Equivalence of Definitions

These and several other definitions turn out to be equivalent under suitable conditions. This was recognized by Rinehart, who wrote in 1955

“There have been proposed in the literature since 1880 eight distinct definitions of a matric function, by Weyr, Sylvester and Buchheim, Giorgi, Cartan, Fantappiè, Cipolla, Schwerdtfeger and Richter … All of the definitions except those of Weyr and Cipolla are essentially equivalent.”

## Computing a Matrix Function in MATLAB

In MATLAB, matrix functions are distinguished from componentwise array evaluation by the postfix “m” on a function name. Thus expm, logm, and sqrtm are matrix function counterparts of exp, log, and sqrt. Compare the example at the start of this article with

>> A = [0 1; 2 3], X = expm(A)
A =
0     1
2     3
X =
5.2892e+00   8.4033e+00
1.6807e+01   3.0499e+01

The MATLAB function funm calculates $f(A)$ for general matrix functions $f$ (subject to some restrictions).

## References

This is a minimal set of references, which contain further useful references within.

## Randsvd Matrices with Large Growth Factors

Sixty years ago James Wilkinson published his backward error analysis of Gaussian elimination for solving a linear system $Ax = b$, where $A$ is a nonsingular $n\times n$ matrix. He showed that in floating-point arithmetic the computed solution $\widehat{x}$ satisfies

$(A+\Delta A) \widehat{x} = b, \qquad \|\Delta A\|_{\infty} \le p(n) \rho_n u \|A\|_{\infty},$

where $u$ is the unit roundoff and $p$ is a low degree polynomial. The term $\rho_n$ is the growth factor, defined by

$\rho_n = \displaystyle\frac{\max_{i,j,k} |a_{ij}^{(k)}|} {\max_{i,j}|a_{ij}|} \ge 1,$

where the $a_{ij}^{(k)}$ are the elements at the $k$th stage of Gaussian elimination. The growth factor measures how much elements grow during the elimination. We would like the product $p(n)\rho_n$ to be of order 1, so that $\Delta A$ is a small relative perturbation of $A$. We therefore need $\rho_n$ not to be too large.

With partial pivoting, in which row interchanges are used to ensure that at each stage the pivot element is the largest in its column, Wilkinson showed that $\rho_n \le 2^{n-1}$ and that equality is possible. Such exponential growth implies a large $\Delta A$ (unless we are lucky), meaning a severe loss of numerical stability. However, seventy years of digital computing experience have shown that $\rho_n$ is usually of modest size in practice. Explaining why this is the case is one of the outstanding problems in numerical analysis.

It is easy to experiment with growth factors in MATLAB. I will use the function

function g = gf(A)
%GF     Approximate growth factor.
%   g = GF(A) is an approximation to the
%   growth factor for LU factorization
%   with partial pivoting.
[~,U] = lu(A);
g = max(abs(U),[],'all')/max(abs(A),[],'all');

It computes a lower bound on the growth factor (since it only considers $k=n$ in the numerator in the definition), but it is entirely adequate for our purposes here. Let’s compute the growth factor for a random matrix of order 10,000 with elements from the standard normal distribution (mean 0, variance 1):

>> rng(1); n = 10000; gf(randn(n))
ans =
6.1335e+01

Growth of 61 is unremarkable for a matrix of this size. Now we try a matrix of the same size generated by the gallery('randsvd') function:

>> A = gallery('randsvd',n,1e6,2,[],[],1);
>> gf(A)
ans =
9.7544e+02

This function generates an $n\times n$ matrix with known singular value distribution and with singular vector matrices that are random orthogonal matrices from the Haar distribution. The parameter 1e6 specifies the 2-norm condition number, while the 2 (the mode parameter) specifies that there is only one small singular value, so the singular values are 1 repeated $n-1$ times and 1e-6. Growth of 975 is exceptional! These matrices have been in MATLAB since the 1990s, but this large growth property has apparently not been noticed before.

It turns out that mode 2 randsvd matrices generate with high probability growth factors of size at least $n/(4 \log n)$ for any condition number and for any pivoting strategy, not just partial pivoting. One way to check this is to randomly permute the columns of $A$ before doing the LU factorization with partial pivoting:

>> gf(A(:,randperm(n)))
ans =
7.8395e+02

Here is a plot showing the maximum over 12 randsvd matrices for each $n$ of the growth factors for three different pivoting strategies, along with the maximum growth factors for partial pivoting for rand and randn matrices. The black curve is $n/(4 \log n)$. This plot emphasizes the unusually large growth for mode 2 randsvd matrices.

What is the explanation for this large growth? It stems from three facts.

• Haar distributed orthogonal matrices have the property that that their elements are fairly small with high probability, as shown by Jiang in 2005.
• If the largest entries in magnitude of $A$ and $A^{-1}$ are both small, in the sense that their product is $\theta \ll 1$, then $A$ will produce a growth factor of at least $1/\theta$ for any pivoting strategy. This was proved by Des Higham and I in the paper Large Growth Factors in Gaussian Elimination with Pivoting.
• If $W$ is an orthogonal matrix generating large growth then a rank-1 perturbation of 2-norm at most 1 tends to preserve the large growth.

For full details see the new EPrint Random Matrices Generating Large Growth in LU Factorization with Pivoting by Des Higham, Srikara Pranesh and me.

Is growth of order $n$ a problem in practice? It can be for two reasons.

• The largest dense linear systems $Ax = b$ solved today are of dimension $n = 10^7$. If we work in single precision then $nu \approx 1$ and so LU factorization can potentially be completely unstable if there is growth of order $n$.
• For IEEE half precision arithmetic growth of order $n$ will cause overflow once $n$ exceeds $10^5 / \max_{i,j} |a_{ij}|$. It was overflow in half precision LU factorization on randsvd matrices that alerted us to the large growth.
Posted in research | Tagged , | 2 Comments

## What Is Bfloat16 Arithmetic?

Bfloat16 is a floating-point number format proposed by Google. The name stands for “Brain Floating Point Format” and it originates from the Google Brain artificial intelligence research group at Google.

Bfloat16 is a 16-bit, base 2 storage format that allocates 8 bits for the significand and 8 bits for the exponent. It contrasts with the IEEE fp16 (half precision) format, which allocates 11 bits for the significand but only 5 bits for the exponent. In both cases the implicit leading bit of the significand is not stored, hence the “+1” in this diagram:

The motivation for bfloat16, with its large exponent range, was that “neural networks are far more sensitive to the size of the exponent than that of the mantissa” (Wang and Kanwar, 2019).

Bfloat16 uses the same number of bits for the exponent as the IEEE fp32 (single precision) format. This makes conversion between fp32 and bfloat16 easy (the exponent is kept unchanged and the significand is rounded or truncated from 24 bits to 8) and the possibility of overflow in the conversion is largely avoided. Overflow can still happen, though (depending on the rounding mode): the significand of fp32 is longer, so the largest fp32 number exceeds the largest bfloat16 number, as can be seen in the following table. Here, the precision of the arithmetic is measured by the unit roundoff, which is $2^{-t}$, where $t$ is the number of bits in the significand.

Note that although the table shows the minimum positive subnormal number for bfloat16, current implementations of bfloat16 do not appear to support subnormal numbers (this is not always clear from the documentation).

As the unit roundoff values in the table show, bfloat16 numbers have the equivalent of about three decimal digits of precision, which is very low compared with the eight and sixteen digits, respectively, of fp32 and fp64 (double precision).

The next table gives the number of numbers in the bfloat16, fp16, and fp32 systems. It shows that the bfloat16 number system is very small compared with fp32, containing only about 65,000 numbers.

The spacing of the bfloat16 numbers is large far from 1. For example, 65280, 65536, and 66048 are three consecutive bfloat16 numbers.

At the time of writing, bfloat16 is available, or announced, on four platforms or architectures.

• The Google Tensor Processing Units (TPUs, versions 2 and 3) use bfloat16 within the matrix multiplication units. In version 3 of the TPU the matrix multiplication units carry out the multiplication of 128-by-128 matrices.
• The NVIDIA A100 GPU, based on the NVIDIA Ampere architecture, supports bfloat16 in its tensor cores through block fused multiply-adds (FMAs) $C + A*B$ with 8-by-8 $A$ and 8-by-4 $B$.
• Intel has published a specification for bfloat16 and how it intends to implement it in hardware. The specification includes an FMA unit that takes as input two bfloat16 numbers $a$ and $b$ and an fp32 number $c$ and computes $c + a*b$ at fp32 precision, returning an fp32 number.
• The Arm A64 instruction set supports bfloat16. In particular, it includes a block FMA $C + A*B$ with 2-by-4 $A$ and 4-by-2 $B$.

The pros and cons of bfloat16 arithmetic versus IEEE fp16 arithmetic are

• bfloat16 has about one less (roughly three versus four) digit of equivalent decimal precision than fp16,
• bfloat16 has a much wider range than fp16, and
• current bfloat16 implementations do not support subnormal numbers, while fp16 does.

If you wish to experiment with bfloat16 but do not have access to hardware that supports it you will need to simulate it. In MATLAB this can be done with the chop function written by me and Srikara Pranesh.

## References

This is a minimal set of references, which contain further useful references within.

## What Is the Matrix Exponential?

The exponential of a square matrix $A$ is defined by the power series (introduced by Laguerre in 1867)

$\mathrm{e}^A = I + A + \displaystyle\frac{A^2}{2!} + \frac{A^3}{3!} + \cdots.$

That the series converges follows from the convergence of the series for scalars. Various other formulas are available, such as

$\mathrm{e}^A = \displaystyle\lim_{s\to \infty} (I+A/s)^s.$

The matrix exponential is always nonsingular and $(\mathrm{e}^A)^{-1} = \mathrm{e}^{-A}$.

Much interest lies in the connection between $\mathrm{e}^{A+B}$ and $\mathrm{e}^A \mathrm{e}^B$. It is easy to show that $\mathrm{e}^{A+B} = \mathrm{e}^A \mathrm{e}^B$ if $A$ and $B$ commute, but commutativity is not necessary for the equality to hold. Series expansions are available that relate $\mathrm{e}^{A+B}$ to $\mathrm{e}^A \mathrm{e}^B$ for general $A$ and $B$, including the Baker–Campbell–Hausdorff formula and the Zassenhaus formula, both of which involve the commutator $[A,B] = AB - BA$. For Hermitian $A$ and $B$ the inequality $\mathrm{trace}(\mathrm{e}^{A+B}) \le \mathrm{trace}(\mathrm{e}^A \mathrm{e}^B)$ was proved independently by Golden and Thompson in 1965.

Especially important is the relation

$\mathrm{e}^A = \bigl(\mathrm{e}^{A/2^s}\bigr)^{2^s},$

for integer $s$, which is used in the scaling and squaring method for computing the matrix exponential.

Another important property of the matrix exponential is that it maps skew-symmetric matrices to orthogonal ones. Indeed if $A = - A^T$ then

$\bigl(\mathrm{e}^A\bigr)^{-1} = \mathrm{e}^{-A} = \mathrm{e}^{A^T} = \bigl(\mathrm{e}^A\bigr)^T.$

This is a special case of the fact that the exponential maps elements of a Lie algebra into the corresponding Lie group.

The matrix exponential plays a fundamental role in linear ordinary differential equations (ODEs). The vector ODE

$\displaystyle\frac{dy}{dt} = A y, \quad y(0) = c$

has solution $y(t) = \mathrm{e}^{At} c$, while the solution of the ODE in $n \times n$ matrices

$\displaystyle\frac{dY}{dt} = AY + YB, \quad Y(0) = C$

is $Y(t) = \mathrm{e}^{At}Ce^{Bt}$.

In control theory, the matrix exponential is used in converting from continuous time dynamical systems to discrete time ones. Another application of the matrix exponential is in centrality measures for nodes in networks.

Many methods have been proposed for computing the matrix exponential. See the references for details.

## References

This is a minimal set of references, which contain further useful references within.

## Related Blog Posts

Posted in what-is | Tagged | 1 Comment

## Simulating Low Precision Floating-Point Arithmetics in MATLAB

At least five floating-point arithmetics are available in mainstream hardware:

• the IEEE double precision (fp64), single precision (fp32), and half precision (fp16) formats,
• bfloat16, and
• tf32, introduced in the recently announced NVIDIA A100, which uses the NVIDIA Ampere GPU architecture.

Only fp32 and fp64 are available on current Intel processors and most programming environments support only these two precisions. For many of us it is therefore necessary to simulate other precisions if we wish to explore their properties and the behavior of algorithms when run at these precisions.

Srikara Pranesh and I have written a MATLAB function chop that allows arithmetic of arbitrary precision and exponent range (both lower than fp64) to be simulated. The idea is very simple: variables are stored in fp32 or fp64 and then rounded to the required precision and exponent range after each operation by calling chop. The precision and range, as well as various other parameters, can be changed at any point in the computation.

Our approach can be compared with defining a new data type and storage format and converting to and from single precision or double precision in order to carry out the computations. Our approach avoids the overheads of repeatedly converting to and from the special storage format. It gains speed and generality at the cost of requiring chop calls to be inserted in the code.

Chop is available from this GitHub repository and from MathWorks File Exchange. The motivation behind chop is described in this open access paper in SIAM. J. Sci. Comput. We have made use of chop to simulate fp16 and bfloat16 in several recent papers, including on mixed precision algorithms, probabilistic rounding error analysis, and fast and accurate summation.

## How to Use Chop

A call to chop has the form chop(X,options), where the structure options spceifies how to round the elements of the real array X to a lower precision arithmetic. X should be a single precision or double precision array and the output will have the same type. The structure options controls various aspects of the rounding, as we now explain.

1. The arithmetic format is specified by options.format, which is one of
• ‘b’, ‘bfloat16’: bfloat16,
• ‘h’, ‘half’, ‘fp16’: IEEE half precision (the default),
• ‘s’, ‘single’, ‘fp32’: IEEE single precision,
• ‘d’, ‘double’, ‘fp64’: IEEE double precision,
• ‘c’, ‘custom’: custom format.

In the last case the (base 2) format is defined by options.params, which is a 2-vector [t,emax], where t is the number of bits in the significand (including the hidden bit) and emax is the maximum value of the exponent.

2. options.subnormal specifies whether subnormal numbers are supported (if they are not, subnormals are flushed to zero):
• 0 = do not support subnormals (the default for bfloat16),
• 1 = support subnormals (the default for fp16, fp32, and fp64).
3. The form of rounding is specified by options.round:
• 1: round to nearest using round to even last bit to break ties (the default);
• 2: round towards plus infinity (round up);
• 3: round towards minus infinity (round down);
• 4: round towards zero;
• 5: stochastic rounding—round to the next larger or next smaller floating-point number with probability proportional to 1 minus the distance to those floating-point numbers;
• 6: stochastic rounding—round to the next larger or next smaller floating-point number with equal probability.

For stochastic rounding, exact floating-point numbers are not changed.

4. If options.flip = 1 (default 0) then each element of the rounded result has, with probability options.p (default 0.5), a randomly chosen bit in its significand flipped. This option is useful for simulating soft errors.
5. If options.explim = 0 (default 1) then emax (the maximal exponent) for the specified arithmetic is ignored, so overflow, underflow, or subnormal numbers will be produced only if necessary for the data type of X. This option is useful for exploring low precisions independent of range limitations.

To avoid repeatedly passing the options argument, one can issue a call chop([],options) and the subsequent calls will use the previously given options parameters. The value of options can be inspected with [~,options] = chop.

In the rest of this post we give examples of the use of chop.

## Fp16 versus Bfloat

The code

x = 1/3;
opt_h.format = 'h'; xh = chop(x,opt_h); err_h = (x - xh)/x
opt_b.format = 'b'; xb = chop(x,opt_b); err_b = (x - xb)/x
x = 70000;
xh = chop(x,opt_h)
xb = chop(x,opt_b)

produces the output

err_h =
2.4414e-04
err_b =
-1.9531e-03
xh =
Inf
xb =
70144

which illustrates the higher precison but smaller range of fp16 compared with bfloat16.

## Rounding Modes

The code

opt.format = 'h'; x = 0.1; rng(2)
for k = 1:6
opt.round = k;
y(k) = chop(x,opt);
end
errors = y - x
diff_y = diff(y)

rounds 0.1 to fp16 using six different rounding modes. The output is

errors =
-2.4414e-05   3.6621e-05  -2.4414e-05  -2.4414e-05  -2.4414e-05
3.6621e-05
diff_y =
6.1035e-05  -6.1035e-05            0            0   6.1035e-05

Rounding maps a number to the next larger or next smaller floating-point number, so there are only two possibilties for the error, and they have opposite signs. The diff_y values are consistent with the spacing of the fp16 floating-point numbers around 0.1, which is $2^{-11}/8$.

## Matrix Multiplication

The overheads of chop can be minimized by choosing a suitable implementation of a computation. Matrix multiplication provides a good example. The multiplication C = A*B of $n$-by-$n$ matrices involves $2n^3$ floating-point operations, each of which needs to be chopped. The following code uses only $2n$ calls to chop, because it processes an outer product in a single call, taking advantage of the fact that chop can take a matrix argument.

A = rand(m,n); B = rand(n,p);
C = zeros(m,p);
for i = 1:n
C = chop(C + chop(A(:,i)*B(i,:)));
end

## Tf32 Arithmetic

The recently announced NVIDIA A100 GPU supports a new floating-point format called tf32 in its tensor cores. Tf32 is a 19-bit format that has 11 bits in the significand (including the hidden bit) and 8 bits in the exponent. So it has the precision of fp16 and the range of bfloat16, effectively combining the best features of these two formats. We can simulate tf32 using the custom format of chop. The code

opt.format = 'custom';
% tf32. Significand: 10 bits plus 1 hidden, exponent: 8 bits.
t = 11; emax = 127;
opt.params = [t,emax];
x = 1/3;
xt = chop(x,opt); err_t = (x - xt)/x
x = 70000;
xt = chop(x,opt)

produces the output

err_t =
2.4414e-04
xt =
70016

The error is the same as for fp16 but the number 70000 now rounds to a finite floating-point number.


text:
I edited an image with a curves adjustment on all the channels (RGB, red, green and blue). Now I want to know how those curves looked like, so I can use them again/edit them. But I made this edit a long time ago so I dont remember.

I have tried simply copying the original edit and managed to get it pretty similar. But not similar enough. The curves were very complex and that together with that I used all 4 channels, makes it very hard to get right.

I thought I had cracked it when I tried to split up the image in its red, green and blue channels. My idea was that it would be easier to copy if I focused on doing one channel at a time in black and white. Then I could save them as preset in Photoshop, and apply them combined to the original. But Ive just spent more than 2 hours getting not close enough on the red channel.

Im not very knowledgeable with curves. Do anyone have an idea on how to pull this off in practice? Could it be possible to "reverse engineer" it somehow?

• If you do not have the original image, and the processed one it is impossible. Do you have both? - And what do you mean you processed 4 channels if you are working in RGB? – Rafael Nov 25 '16 at 17:09
• I have both the original and the edit. Red channel + green channel + blue channel + the master RBG channel. Thats 4. – ismathme Nov 26 '16 at 17:58

If you'd have used a single curve per channel, and you have your input and output images, you could derive the curves back using some image processing.

But since you used both the RGB curve and the individual R, G and B curves, that is (I think) impossible (or exceedingly hard), since you won't be able to tell how much the RGB curve or the individual curves each contributed to the output.

• Since he has both the original and the adjusted one then no this is not at all hard just plot the value difference in a graph and your done with it – joojaa Mar 6 '18 at 12:15
• Since the poster mentioned they adjusted both the RGB curve and the R/G/B curves, it's harder (but assume not impossible?) to get back the original curves (since each channel will have been affected by both the RGB curve and the individual channel curve), but I assume you could get back something useful from the plot, yes. – Dan Burzo Mar 6 '18 at 12:47
• Well it really depends on what tools you have at your disposal i tried to do this in mathematica and i have allready a solution that works fine – joojaa Mar 6 '18 at 12:51
• added an answer on the subject while i was testing this. – joojaa Mar 6 '18 at 13:55
• Cool! Does it work with changing both the RGB channel and the individual R/G/B channels? That's the part I was trying to clarify: I can't intuit how'd you extract the original curves if there are two curves per channel that affect the output. – Dan Burzo Mar 6 '18 at 14:02

Yes you can do this. However it really depends on many factors. First of you can just simply use some scientific plotting tool to plot each value in 2D where x is original value and y is new pictures value.

I am using mathematica for this, quick and dirty mathematica code looks as follows (but you could choose any other tool):

src = ColorSeparate[Import["c:\\temp\\test.png"]];
r = Flatten[ImageData [src[[1]]]];
rtp = Transpose[{r, r2}];
g = Flatten[ImageData [src[[2]]]];
gtp = Transpose[{g, g2}];
b = Flatten[ImageData [src[[3]]]];
btp = Transpose[{b, b2}];
Show[ListPlot[rtp, AspectRatio -> 1, PlotStyle -> Red],
ListPlot[gtp, AspectRatio -> 1, PlotStyle -> Green],
ListPlot[btp, AspectRatio -> 1, PlotStyle -> Blue]]

So for test ive used following curves as the adjusted image (image used):

And mathematica gives me:

Which is pretty close, or possibly exactly the same value. Making a curves adjustment file out of this wouldn't be too hard*. However a word of warning! If your images are lossy, like jpeg compressed, then this will not work. Example with jpegs:

Now the plot points are all over the place and youd have to try to filter the data somehow. So the take here is that if the image has something other than just a curves adjustment then the thing just blows up.

* If the data was read in with the code above then the curves export is:

file = FileNameJoin[{ "C:\\temp", "out.amp"}];
BinaryWrite[
file, (SortBy[DeleteDuplicates[rtp], First] // Transpose //
Last) *255 // IntegerPart];
BinaryWrite[
file, (SortBy[DeleteDuplicates[gtp], First] // Transpose //
Last) *255 // IntegerPart];
BinaryWrite[
file, (SortBy[DeleteDuplicates[btp], First] // Transpose //
Last) *255 // IntegerPart];
Close[file];

the amp file can then be read to a photoshop curves layer for example. And it produces same image in my test case. The code is not really performace oriented but it works.

• @dan Hey, thanks guys! But to be honest i dont even remember what image it was i wanted to do this with haha. I read your answers and it was interesting though. I dont really have a way to test it, but should i mark joojaa's answer as a solution anyways? – ismathme Mar 9 '18 at 13:33
• @ismathme its your question. BUt out of curiosty what tools do you have. Theres no reason why this couldnt be ported to say python and numpy/scipy for example – joojaa Mar 9 '18 at 20:16
• Well i dont really have any tools other than photoshop. I do have a little programming knowledge, so i was able to follow your code in broad strokes. But no, i dont really have a way to actually execute this lol – ismathme Mar 10 '18 at 13:46


text:
The logo of the show “The Fresh Prince of Bel-Air”, originally aired on NBC. Copyright Warner Bros.

# The Pharo Prince of Bel-Air

Recently I have been aiming to explore the Pharo programming language. The close combination between the language and the Integrated Development Environment (IDE), makes it quite unlike most other languages that I have encountered. You can, for example, inspect every program element in a visual environment, even on a currently running program. While exploring, I spent much time using the Playground feature, where you easily start executing, and interacting with, small snippets of code.

While exploring the language, the intro-song to the tv-show “The Fresh Prince of Bel-Air” got stuck in my head. The song tells the story of the title character, Will Smith, on how he ended up being “The Fresh Prince of Bel-Air”. To paraphrase the song, Will was born and raised in Philadelphia, spending much of his time on the playground. On one occasion, while playing basketball, he got into a fight with some troublemakers. This experience got his mom scared so he was sent off to live with his aunt and uncle in Bel Air. After a flight and a taxi ride, he arrives in Bel Air at his uncles and aunts house, to become “the Prince of Bel-Air”.

In this article, let’s explore Pharo together with the story from “The Fresh Prince of Bel Air”. If one wants to follow along, feel free to download and setup Pharo. The language has a really smooth process for getting started, as one can download and install the Pharo Launcher and be pretty much set.

It is difficult to do the language justice by text, because it is so intertwined with the IDE and image based development. I will try to illustrate things as best as I can.

Once we have installed and started the Pharo Launcher the first step is to create a Pharo image for our project. One can think of a Pharo image, as a store of all the objects created by the program. This not only includes what would normally be represented as files of source code of other languages, but also the objects of the running program, as well as the state of the editor. As an example, suppose we create an application for taking notes, and while developing, we create two notes to experiment and play around with. We can now save the Pharo image including these objects. The next time we start the image our notes will be there, alongside with the state of the IDE as we have left it. This can help tremendously when developing or debugging a new application.

To come back to our “The Fresh Prince of Bel-Air” example, first we are going to create an image which will contain our project (see a screenshot of this below).

In the Pharo launcher, we can click on the “New” button and get a dialog window for the image that we want to create. We base our image on the last stable Pharo version, which is at the time of this writing version 8.0, enter a name and a short description for our image and click “Create Image”.

In a short order our image is then initialized, and we are greeted in our IDE by a welcome message, and some guides and documentation in a menu. I highly recommend the “ProStef” tutorial from this set of guides.

Our next step is the playground which, as mentioned at the start of the article, is a place for experimentation and playing around with code that we are going to develop. One can open this Playground from the menu “Tools” on the top of the screen by selecting “Playground” or by the shortcut CTRL+O+W.

Note that the shortcuts might be different depending on settings and operating systems, as the command key might be either CTRL, CMD or ALT. Thankfully the IDE shows shortcuts very neatly next to the clickable menu options, so they should be easy to figure out and learn.

Once the Playground is open we have a canvas where we can start to write and execute our code.

Before I get into the actual creation of the example I will briefly try to summarize how Pharo as a language works. In Pharo everything is an object and any computation is done by sending messages to objects. Creating such objects is done through specific objects known as classes, which define structure and behaviour of the objects they create.

These classes are organized into hierarchies. Hierarchies allow you to reuse behaviour. A class lower in the hierarchy, the so-called subclass can inherit the behaviour of a class higher in the hierarchy, called a superclass. At the root of this hierarchy is the class named Object, which ensures as we noted that anything we create is by definition and behaviour an Object.

We call the objects created by the class as their instances. The methods of an object specify how the messages an object receives should be handled. When an object receives a message it aims to find the method that has the same name as the message to handle it. An object also encapsulates data, such as dates, numbers, strings and other objects, based on what the object represents into variables. It can use this data to help handle the messages it receives.

An object can refer to itself by using the special identifier self. In order to create a new object of a class, the initialize message needs to be sent which gets handled by the initialize method, though often the message new is used as a shorthand for this.

For people with some experience with another Object Oriented programming language the above description should sound familiar. One thing that stands out for me on a first glance is how consistently the use of objects and message passing is applied within Pharo. Everything is an object and interaction seems to always happen by message passing. The above summary should hopefully be enough to follow the examples below but for a more in depth look I highly recommend the Learning Object-Oriented Programming, Design with TDD in Pharo book.

Now given the above description of Pharo, we want to apply it by creating a representation of what happens in the intro song in “The Fresh Prince of Bel-Air” by sending messages. In particular we want to represent the character of “Will Smith” in the intro as an object and the changes in location and activity based on the intro song.

So let’s start off by creating a new Object representing Will. We aim to do this by instantiating a class called FPWillSmith by typing the following into the Playground.

will := FPWillSmith new.

Here we want to send the New message to the class FPWillSmith to create a new object of this class and assign it to the variable will. Let us quickly go over the syntax elements for this statement.

The operator := is used in Pharo for assignment, and in this case to assign our new object instance to the will variable. We also send to class FPWillSmith, functioning as the receiver, the message new to create a new instance of it. Messages like this that are sent to an object without any additional information are called unary messages. These are always in the form of receiver message. Lastly, statements in Pharo are separated/ended by a .. This can be omitted in many cases but we keep it around in our examples to show its use.

To execute this statement we can right click it and select ‘Do It’ or use the short-cut for this CTRL+D.

After doing this we will get the following dialog:

Well, whoops it seems we kind of forgot to define our class FPWillSmith, so we get a notification for it. The reason why I am highlighting this error is to show that unlike many other languages, where you would just get an error at either compile or runtime, the system brings up a dialog to help us solve this issue. The only other language where I experienced a similar dialog was with CommonLisp.

So lets click “Define a new class” and give a class definition as follows:

A few things to note for this code in this picture:

Object subclass: #FPWillSmith
instanceVariableNames: 'name location activity'
classVariableNames: ''
category: FreshPrinceDemo

Much like every other part of Pharo, the creation of a new class also happens through message passing. Here a keyword message is sent to the Object class to define our new class. Unlike unary messages, these messages use one or more keywords to send additional object information along with the message. These are defined in the form of receiver keyword1: object1 keyword2: object2 ... where the keyword and object pairs make up the message to the receiving object.

We will go over the objects sent in the message one by one. The subclass part of the message defines the name of our new class FPWillSmith. The # in front of FPWillSmith indicates a symbol that is globally unique, unlike a string, which is part of the requirement for this subclass message. Next is defining the name of the instance variables using the keyword instanceVariablesNames, that will represent the elements of an instance of FPWillSmith that we would like to use. There are also classVariableNames that describe such elements for the class itself as opposed to the instances, but here we do not define any. Finally, there is a category defined as FreshPrinceDemo, that is essentially the name of the package that we are currently defining.

Now that we have created our FPWillSmith class, we can now take our original statement:

will := FPWillSmith new.

and execute it, by right clicking it and selecting ‘Do It’, or using a shortcut such as CTRL+D. As you can see, we now don’t receive errors.

What we would like to do is to double check, whether the object we created is what we want. Pharo has a really neat feature for inspecting any object. We can right-click this statement in the playground and select ‘Inspect It’ from the menu, or using the shortcut CTRL+I for giving us the object inspector. This should give us a window such as this:

As we can see from the inspector, while the object has all the instance variables that we wanted (name, location, activity), the value for this is nil which is not what we want.

In order to make sure there are meaningful values for these variables for a new object, we need to define our own initialize method that will get called for every new object.

First we are going to look up our class in the code browser. This can be done by clicking the menu options Tools and then Systems Browser, or by using the short-cut CTRL+O+B.

The system browser can be used to find, organize and edit classes and methods for all packages within the Pharo image. It is divided up into four columns. The leftmost column lists the packages. We can use this to filter for our new package named ‘FreshPrinceDemo’. The second column on the left details the classes in the selected package. For our new package we should only have one class here: FPWillSmith. The third column lists all the protocols of the selected class. Protocols are a way to group together related methods, which makes for easier browsing and search. Finally the fourth column from the left lists all the methods, for the given class or selected protocol.

The panel below these four columns is also very interesting. If one selects a package, it shows in a tab a template for defining a class. If one selects a class, it shows a tab for a template on how to define new methods. It also has a tab for the class definition itself .

In order to proceed let’s use the method definition template to define our new constructor, as follows:

initialize
name := 'Will Smith'.

which we can fill into the new instance side method window below, instead of the template code.

What this definition does is to create a constructor, that sets the values for name, location and activity to values that makes sense, given where Will is and what he is doing in the beginning of the “The Fresh Prince of Bel-Air” intro song.

Make sure to accept these values, using right click on the panel and selecting Accept or by pressing the shortcut CTRL+S (you might also need to save your name as the author for this piece of code in a dialog along the way).

We will now have the initialization method defined:

We can check whether our new constructor has worked by inspecting the statement we initially created in the Playground:

Now the will representation of the character has the right name, location and activity, but we want to express a bit more of the intro song. We will do this by sending messages to the will object after which we expect that the location and the activity changes accordingly.

Let’s assume we instead have the following in the Playground:

will := FPWillSmith new.
will getsInTrouble.

As we can now expect, executing (‘Doing It’) of the second statement will not go well as we did not yet define a way for the object to handle this message. In fact this message will bring up the debugger because the object does not know how to handle the message getsInTrouble.

We can fix this by defining a method, named getsInTrouble. We can even do this from the debugger by clicking on the Create button at the top of the debugger to define getsInTrouble method for the FPWillSmith class. After getting a popup for picking a protocol, which we do not need to do for this tutorial, we get the window showing a placeholder definition for the method:

getsInTrouble
self shouldBeImplemented.

We want to replace this with the following definition, that updates the activity of the object:

getsInTrouble
activity := 'In trouble'.

We have now created the method from inside the debugger, while the execution of the code was put on hold! With the definition saved we can even click the Proceed button on the top of the debugger to continue from where we halted.

Now the code correctly executes and when we next inspect the last statement have the activity properly updated.

For the final piece we will define a method for sending the object to aunt and uncle, with the following definition:

sentToAuntAndUncle
activity := 'Being the Prince of Bel-Air'.
location := 'Bel-Air'.

If one follows along, we will leave it up to you on where to define this method from, the debugger or the system browser.

Once we have defined it we can put the following in the Playground:

will := FPWillSmith new.
will getsInTrouble.
will sentToAuntAndUncle.

and inspect our object:

That is it for this short tutorial. I hope I have shown a small glimpse of why Pharo is such an interesting language.

We have started with no definitions at all, and slowly built up what we wanted to do in a very organic way. We just started putting the statements down in the Playground for passing messages to objects, and created the definitions as we went along. Whenever something has given an error in the debugger or did not show up to our liking in the inspector, we slowly added more code. In the process we have evolved an object that initially represented no aspect of Will’s character, to something that had a name, location and activity, including; “Playing basketball” in a “West Philadelphia playground”.

And at the end we are now finally there, having the object’s activity denoted as Being the Prince of Bel-Air.


text:
# Special Triangles Question

• June 1st 2011, 05:51 AM
starshine84
Special Triangles Question
I am looking for help with this question from my grade 11 advanced functions course:

csc θ = -17/15 where 270° ≤ θ ≤ 360° and
cot β = - 3/4 where 90° ≤ β ≤ 180°.

Find the exact value of sin(θ + β). Show all the work.
• June 1st 2011, 05:54 AM
Prove It
Better yet, why don't you show all of YOUR work, i.e. what you have tried and where exactly you are stuck... Then you might get the help you need...
• June 1st 2011, 06:03 AM
starshine84
This what I have done:

cscθ = -17/15 where 270° ≤ θ ≤ 360°

sinθ = 1/cscθ
sinθ = 1/(-17/15)
sinθ = 1/1 / (-17/15)
sinθ = 1/1 * (-15/17)
sinθ = -15/17
θ = sin^-1 (-15/17)

** I don't know what to do from here because if I solve for θ completely I won't have an EXACT value **

cot β = - 3/4 where 90° ≤ β ≤ 180°

tanβ = 1/cotβ
tanβ = 1/(-3/4)
tanβ = 1/1 / (-3/4)
tanβ = 1/1 * (-4/3)
tanβ = -4/3
β = tan^-1 (-4/3)

** I don't know what to do from here because if I solve for β completely I won't have an EXACT value **

From this point I'm not sure how to find an exact value for sin(θ + β) because I do not have values for θ or β
• June 1st 2011, 06:19 AM
Prove It
For starters, you should know that $\displaystyle \sin{(\theta \pm \beta)} \equiv \sin{(\theta)}\cos{(\beta)} \pm \cos{(\theta)}\sin{(\beta)}$.

You don't need to solve for $\displaystyle \theta$ or $\displaystyle \beta$, you only need to solve for $\displaystyle \sin{(\theta)}, \cos{(\theta)}, \sin{(\beta)}$ and $\displaystyle \cos{(\beta)}$. You will need to make use of the Pythagorean Identity.
• June 1st 2011, 07:46 AM
starshine84
I knew about the sum and difference formulas and that I would have to use the sum formula for sin in order to solve the final step, but I'm not sure what you mean by pythagorean identities. I'm not sure how to solve for sin(\theta), cos(\theta), sin(\beta), or cos(\beta).
• June 1st 2011, 08:12 AM
Plato
Quote:

Originally Posted by starshine84
This what I have done:
cscθ = -17/15 where 270° ≤ θ ≤ 360°
cot β = - 3/4 where 90° ≤ β ≤ 180°
** I don't know what to do from here because if I solve for β completely I won't have an EXACT value

From that given you know that $\theta\in IV$.
So $\sin \left( \theta \right) = \frac{{ - 15}}{{17}}\;\& \,\cos \left( \theta \right) = \sqrt {1 - \left( {\frac{15}{17}} \right)^2 }=\frac{8}{{17}}$

Moreover, $\beta\in II$ so $\sin \left( \beta \right) = \frac{4}{5}\;\& \,\cos \left( \beta \right) = \frac{{ - 3}}{5}$.
• June 1st 2011, 08:41 AM
starshine84
So would it then be...

sin(\theta+\beta) = sin\theta cos\beta + cos\theta sin\beta
sin(\theta+\beta) = (-15/17)(-3/4) + (8/17)(4/5)
sin(\theta+\beta) = (45/68) + (32/85)
sin(\theta+\beta) = 225/340 + 128/340
sin(\theta+\beta) = 353/340

??
• June 1st 2011, 02:58 PM
Prove It
Looks good to me :)
• June 2nd 2011, 04:05 AM
starshine84
That's awesome, thank you so much for all of your help!