text:
Warning

This is A PREVIEW for NEST 3.0 and NOT an OFFICIAL RELEASE! Some functionality may not be available and information may be incomplete!

# Part 1: Neurons and simple neural networks¶

## Introduction¶

In this section we cover the first steps in using PyNEST to simulate neuronal networks. When you have worked through this material, you will know how to:

• start PyNEST

• create neurons and stimulating/recording devices

• query and set their parameters

• connect them to each other or to devices

• simulate the network

• extract the data from recording devices

More advanced examples can be found at Example Networks, or have a look at at the source directory of your NEST installation in the subdirectory: pynest/examples/.

## PyNEST - an interface to the NEST simulator¶

Figure 1 Python Interface Figure. The Python interpreter imports NEST as a module and dynamically loads the NEST simulator kernel (pynestkernel.so). The core functionality is defined in hl_api.py. A simulation script of the user (mysimulation.py) uses functions defined in this high-level API. These functions generate code in SLI (Simulation Language Interpreter), the native language of the interpreter of NEST. This interpreter, in turn, controls the NEST simulation kernel.

The NEural Simulation Tool (NEST: www.nest-initiative.org) 1 is designed for the simulation of large heterogeneous networks of point neurons. It is open source software released under the GPL licence. The simulator comes with an interface to Python 2. Figure 1 illustrates the interaction between the user’s simulation script (mysimulation.py) and the NEST simulator. Eppler et al. 3 contains a technically detailed description of the implementation of this interface and parts of this text are based on this reference. The simulation kernel is written in C++ to obtain the highest possible performance for the simulation.

You can use PyNEST interactively from the Python prompt or from within ipython. This is very helpful when you are exploring PyNEST, trying to learn a new functionality or debugging a routine. Once out of the exploratory mode, you will find it saves a lot of time to write your simulations in text files. These can in turn be run from the command line or from the Python or ipython prompt.

Whether working interactively, semi-interactively, or purely executing scripts, the first thing that needs to happen is importing NEST’s functionality into the Python interpreter.

import nest

It should be noted, however, that certain external packages must be imported before importing nest. These include scikit-learn and SciPy.

from sklearn.svm import LinearSVC
from scipy.special import erf

import nest

As with every other module for Python, the available functions can be prompted for.

dir(nest)

One such command is nest.Models() or in ipython nest.Models?, which will return a list of all the available models you can use. If you want to obtain more information about a particular command, you may use Python’s standard help system.

This will return the help text (docstring) explaining the use of this particular function. There is a help system within NEST as well. You can open the help pages in a browser using nest.helpdesk() and you can get the help page for a particular object using nest.help(object).

## Creating Nodes¶

A neural network in NEST consists of two basic element types: nodes and connections. Nodes are either neurons, devices or sub-networks. Devices are used to stimulate neurons or to record from them. Nodes can be arranged with spatial structure to build networks distributed in space - we will get to this later in the course. For now we will work with the default network structure of NEST.

New nodes are created with the command Create, which takes as arguments the model name of the desired node type, and optionally the number of nodes to be created and the initialising parameters. The function returns a NodeCollection of handles to the new nodes, which you can assign to a variable for later use. A NodeCollection is a compact representation of the node handles, which are integer numbers, called ids. Many PyNEST functions expect or return a NodeCollectoin (see command overview). Thus, it is easy to apply functions to large sets of nodes with a single function call.

After having imported NEST and Matplotlib 4, which we will use to display the results, we can start creating nodes. As a first example, we will create a neuron of type iaf_psc_alpha. This neuron is an integrate-and-fire neuron with alpha-shaped postsynaptic currents. The function returns a NodeCollection of the ids of all the created neurons, in this case only one, which we store in a variable called neuron.

import matplotlib.pyplot as plt
import nest
neuron = nest.Create("iaf_psc_alpha")

We can now use the NodeCollection to access the properties of this neuron. Properties of nodes in NEST are generally accessed via Python dictionaries of key-value pairs of the form {key: value}. In order to see which properties a neuron has, you may ask it for its status.

neuron.get()

This will print out the corresponding dictionary in the Python console. Many of these properties are not relevant for the dynamics of the neuron. To find out what the interesting properties are, look at the documentation of the model through the helpdesk. If you already know which properties you are interested in, you can specify a key, or a list of keys, as an optional argument to get:

neuron.get("I_e")
neuron.get(["V_reset", "V_th"])

In the first case we query the value of the constant background current I_e; the result is given as a floating point element. In the second case, we query the values of the reset potential and threshold of the neuron, and receive the result as a dictionary . If get is called on a NodeCollection with more than one element, the returned dictionary will contain lists with the same number of elements as the number of nodes in the NodeCollection. If get is called with a specific key on a NodeCollection with several elements, a list the size of the NodeCollection will be returned.

To modify the properties in the dictionary, we use set. In the following example, the background current is set to 375.0pA, a value causing the neuron to spike periodically.

neuron.set(I_e=375.0)

Note that we can set several properties at the same time by giving multiple comma separated key:value pairs in a dictionary. Also be aware that NEST is type sensitive - if a particular property is of type double, then you do need to explicitly write the decimal point:

neuron.set({"I_e": 375})

will result in an error. This conveniently protects us from making integer division errors, which are hard to catch.

Another way of setting and getting parameters is to ask the NodeCollection directly

neuron.I_e = 376.0
neuron.I_e

Next we create a multimeter, a device we can use to record the membrane voltage of a neuron over time. The property record_from expects a list of the names of the variables we would like to record. The variables exposed to the multimeter vary from model to model. For a specific model, you can check the names of the exposed variables by looking at the neuron’s property recordables.

multimeter = nest.Create("multimeter")
multimeter.set(record_from=["V_m"])

We now create a spikedetector, another device that records the spiking events produced by a neuron.

spikedetector = nest.Create("spike_detector")

A short note on naming: here we have called the neuron neuron, the multimeter multimeter and so on. Of course, you can assign your created nodes to any variable names you like, but the script is easier to read if you choose names that reflect the concepts in your simulation.

## Connecting nodes with default connections¶

Now we know how to create individual nodes, we can start connecting them to form a small network.

nest.Connect(multimeter, neuron)
nest.Connect(neuron, spikedetector)

Figure 2 Membrane potential of integrate-and-fire neuron with constant input current.

Figure 3 Spikes of the neuron.

The order in which the arguments to Connect are specified reflects the flow of events: if the neuron spikes, it sends an event to the spike detector. Conversely, the multimeter periodically sends requests to the neuron to ask for its membrane potential at that point in time. This can be regarded as a perfect electrode stuck into the neuron.

Now we have connected the network, we can start the simulation. We have to inform the simulation kernel how long the simulation is to run. Here we choose 1000ms.

nest.Simulate(1000.0)

Congratulations, you have just simulated your first network in NEST!

## Extracting and plotting data from devices¶

After the simulation has finished, we can obtain the data recorded by the multimeter.

dmm = multimeter.get()
Vms = dmm["events"]["V_m"]
ts = dmm["events"]["times"]

In the first line, we obtain a dictionary with status parameters for the multimeter. This dictionary contains an entry named events which holds the recorded data. It is itself a dictionary with the entries V_m and times, which we store separately in Vms and ts, in the second and third line, respectively. If you are having trouble imagining dictionaries of dictionaries and what you are extracting from where, try first just printing dmm to the screen to give you a better understanding of its structure, and then in the next step extract the dictionary events, and so on.

Now we are ready to display the data in a figure. To this end, we make use of matplotlib and the pyplot module.

import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(ts, Vms)

The second line opens a figure (with the number 1), and the third line actually produces the plot. You can’t see it yet because we have not used plt.show(). Before we do that, we proceed analogously to obtain and display the spikes from the spike detector.

dSD = spikedetector.get("events")
evs = dSD["senders"]
ts = dSD["times"]
plt.figure(2)
plt.plot(ts, evs, ".")
plt.show()

Here we extract the events more concisely by sending the parameter name to get. This extracts the dictionary element with the key events rather than the whole status dictionary. The output should look like Figure 2 and Figure 3. If you want to execute this as a script, just paste all lines into a text file named, say, one-neuron.py . You can then run it from the command line by prefixing the file name with python, or from the Python or ipython prompt, by prefixing it with run.

It is possible to collect information of multiple neurons on a single multimeter. This does complicate retrieving the information: the data for each of the n neurons will be stored and returned in an interleaved fashion. Luckily Python provides us with a handy array operation to split the data easily: array slicing with a step (sometimes called stride). To explain this you have to adapt the model created in the previous part. Save your code under a new name, in the next section you will also work on this code. Create an extra neuron with the background current given a different value:

neuron2 = nest.Create("iaf_psc_alpha")
neuron2.set({"I_e": 370.0})

now connect this newly created neuron to the multimeter:

nest.Connect(multimeter, neuron2)

Run the simulation and plot the results, they will look incorrect. To fix this you must plot the two neuron traces separately. Replace the code that extracts the events from the multimeter with the following lines.

plt.figure(2)
Vms1 = dmm["events"]["V_m"][::2] # start at index 0: till the end: each second entry
ts1 = dmm["events"]["times"][::2]
plt.plot(ts1, Vms1)
Vms2 = dmm["events"]["V_m"][1::2] # start at index 1: till the end: each second entry
ts2 = dmm["events"]["times"][1::2]
plt.plot(ts2, Vms2)

Additional information can be found at http://docs.scipy.org/doc/numpy-1.10.0/reference/arrays.indexing.html.

## Connecting nodes with specific connections¶

A commonly used model of neural activity is the Poisson process. We now adapt the previous example so that the neuron receives 2 Poisson spike trains, one excitatory and the other inhibitory. Hence, we need a new device, the poisson_generator. After creating the neurons, we create these two generators and set their rates to 80000Hz and 15000Hz, respectively.

noise_ex = nest.Create("poisson_generator")
noise_in = nest.Create("poisson_generator")
noise_ex.set(rate=80000.0)
noise_in.set(rate=15000.0)

Additionally, the constant input current should be set to 0:

neuron.set(I_e=0.0)

Each event of the excitatory generator should produce a postsynaptic current of 1.2pA amplitude, an inhibitory event of -2.0pA. The synaptic weights can be defined in a dictionary, which is passed to the Connect function using the keyword syn_spec (synapse specifications). In general all parameters determining the synapse can be specified in the synapse dictionary, such as "weight", "delay", the synaptic model ("synapse_model") and parameters specific to the synaptic model.

syn_dict_ex = {"weight": 1.2}
syn_dict_in = {"weight": -2.0}
nest.Connect(noise_ex, neuron, syn_spec=syn_dict_ex)
nest.Connect(noise_in, neuron, syn_spec=syn_dict_in)

Figure 4 Membrane potential of integrate-and-fire neuron with Poisson noise as input.

Figure 5 Spikes of the neuron with noise.

The rest of the code remains as before. You should see a membrane potential as in Figure 4 and Figure 5.

In the next part of the introduction (Part 2: Populations of neurons) we will look at more methods for connecting many neurons at once.

## Two connected neurons¶

Figure 6 Postsynaptic potentials in neuron2 evoked by the spikes of neuron1

There is no additional magic involved in connecting neurons. To demonstrate this, we start from our original example of one neuron with a constant input current, and add a second neuron.

import nest
neuron1 = nest.Create("iaf_psc_alpha")
neuron1.set(I_e=376.0)
neuron2 = nest.Create("iaf_psc_alpha")
multimeter = nest.Create("multimeter")
multimeter.set(record_from=["V_m"])

We now connect neuron1 to neuron2, and record the membrane potential from neuron2 so we can observe the postsynaptic potentials caused by the spikes of neuron1.

nest.Connect(neuron1, neuron2, syn_spec = {"weight":20.0})
nest.Connect(multimeter, neuron2)

Here the default delay of 1ms was used. If the delay is specified in addition to the weight, the following shortcut is available:

nest.Connect(neuron1, neuron2, syn_spec={"weight":20.0, "delay":1.0})

If you simulate the network and plot the membrane potential as before, you should then see the postsynaptic potentials of neuron2 evoked by the spikes of neuron1 as in Figure 6.

## Command overview¶

These are the functions we introduced for the examples in this handout; the following sections of this introduction will add more.

See the Getting Help Section

### Nodes¶

• Create(model, n=1, params=None)

Create n instances of type model. Parameters for the new nodes can be given as params (a single dictionary, or a list of dictionaries with size n). If omitted, the model’s defaults are used.

• get(*params, **kwargs)

Return a dictionary with parameter values for the NodeCollection it is called on. If params is a single string, a list of values is returned instead. params may also be a list of strings, in which case the returned dictionary contains lists of requested values.

• set(params=None, **kwargs)

Set the parameters on the NodeCollection it is called on to params, which may be a single dictionary, or a list of dictionaries of the same size as the NodeCollection. If kwargs is given, it has to be names and values of an attribute as keyword=argument pairs. The values can be single values or list of the same size as the NodeCollection.

### Connections¶

This is an abbreviated version of the documentation for the Connect function, please see NEST’s online help for the full version and Connection Management for an introduction and examples.

• Connect(pre, post, conn_spec=None, syn_spec=None, return_synapsecollection=False)

Connect pre neurons to post neurons. Neurons in pre and post are connected using the specified connectivity ("all_to_all" by default) and synapse type ("static_synapse" by default). Details depend on the connectivity rule. pre - presynaptic neurons, given as a NodeCollection of node IDs post - presynaptic neurons, given as a NodeCollection of node IDs conn_spec - name or dictionary specifying connectivity rule, see below syn_spec - name or dictionary specifying synapses, see below.

#### Connectivity¶

Connectivity is either specified as a string containing the name of a connectivity rule (default: "all_to_all") or as a dictionary specifying the rule and rule-specific parameters (e.g. "indegree"), which must be given. In addition switches allowing self-connections ("allow_autapses", default: True) and multiple connections between a pair of neurons ("allow_multapses", default: True) can be contained in the dictionary.

#### Synapse¶

The synapse model and its properties can be inserted either as a string describing one synapse model (synapse models are listed in the synapsedict) or as a dictionary. If no synapse model is specified the default model "static_synapse" will be used. Available keys in the synapse dictionary are "synapse_model", "weight", "delay", "receptor_type" and parameters specific to the chosen synapse model. All parameters are optional and if not specified will use the default values determined by the current synapse model. "synapse_model" determines the synapse type, taken from pre-defined synapse types in NEST or manually specified synapses created via CopyModel(). All other parameters can be scalars or distributions. In the case of scalar parameters, all keys take doubles except for "receptor_type" which has to be initialised with an integer. Distributed parameters are initialised with a Parameter with distribution-specific arguments (such as "mean" and "std").

### Simulation control¶

• Simulate(t)

Simulate the network for t milliseconds.

## References¶

1

Gewaltig MO. and Diesmann M. 2007. NEural Simulation Tool. 2(4):1430.

2

Python Software Foundation. The Python programming language, 2008. http://www.python.org.

3

Eppler JM et al. 2009 PyNEST: A convenient interface to the NEST simulator. 2:12. 10.3389/neuro.11.012.2008.

4

Hunter JD. 2007 Matplotlib: A 2d graphics environment. 9(3):90–95.


text:
DESeq2 paired multifactor test
4
1
Entering edit mode
hwillenbrock ▴ 10
@hwillenbrock-6937
Last seen 6.8 years ago
Denmark

Hi,

I am wondering what the correct model formula and analysis approach is with DESeq2 for my multifactor paired design. I have tried to look at previous answers, e.g. DESEq2 comparison with two conditions and paired experiments and DESeq2 multiple factors nested design but I keep getting warnings from DESeq2, so I must be interpreting them wrong, or my design is slightly different.

What I have is 2 different tissue cell cultures, 5 replicates, that have been treated with a chemical. I have RNAseq count data for both the untreated and treated samples, so each treated vs. untreated is a pair:

sample   tissue culture    treatment
1   tissueX 1 unstimulated
2   tissueX 1 stimulated
3   tissueX 2 unstimulated
4   tissueX 2 stimulated
5   tissueX 3 unstimulated
6   tissueX 3 stimulated
7   tissueX 4 unstimulated
8   tissueX 4 stimulated
9   tissueX 5 unstimulated
10  tissueX 5 stimulated
11  tissueY 6 unstimulated
12  tissueY 6 stimulated
13  tissueY 7 unstimulated
14  tissueY 7 stimulated
15  tissueY 8 unstimulated
16  tissueY 8 stimulated
17  tissueY 9 unstimulated
18  tissueY 9 stimulated
19  tissueY 10 unstimulated
20  tissueY 10 stimulated

The genes that I'm interested in identifying, are those where the stimulation effect differs most between the two tissues. If not taking the paired design into account, I can use the model: ~tissue+treatment+tissue:treatment that is accepted without errors by DESeq2. However, I don't get very significant interactions and I'd like to take advantage of the paired design.

If I just analyze tissueX unstimulated versus stimulated, I can add the "culture" factor to the model and it improves the pvalues considerably: ~subject+treatment

Simply adding the subject term to the model formula as for the simple one tissue paired case, like this: ~subject + tissue+treatment+tissue:treatment gives me an error.

Trying to interpolate a previous answer to my case: ~ subject + treatment + tissue:treatment gives me the same error: "the model matrix is not full rank, so the model cannot be fit as specified."

How do I incorporate the paired design into the model formula and extract the correct significance values?

Hanni

> sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-suse-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
[1] ggplot2_1.0.0             DESeq2_1.6.1
[5] GenomicRanges_1.18.1      GenomeInfoDb_1.2.2
[7] IRanges_2.0.0             S4Vectors_0.4.0
[9] Biobase_2.26.0            BiocGenerics_0.12.0

loaded via a namespace (and not attached):
[1] acepack_1.3-3.3      annotate_1.44.0      AnnotationDbi_1.28.0
[4] base64enc_0.1-2      BatchJobs_1.4        BBmisc_1.7
[7] BiocParallel_1.0.0   brew_1.0-6           checkmate_1.5.0
[10] cluster_1.15.3       codetools_0.2-9      colorspace_1.2-4
[13] DBI_0.3.1            digest_0.6.4         fail_1.2
[16] foreach_1.4.2        foreign_0.8-61       Formula_1.1-2
[19] genefilter_1.48.1    geneplotter_1.44.0   grid_3.1.0
[22] gtable_0.1.2         Hmisc_3.14-5         iterators_1.0.7
[25] labeling_0.3         lattice_0.20-29      latticeExtra_0.6-26
[28] locfit_1.5-9.1       MASS_7.3-31          munsell_0.4.2
[31] nnet_7.3-8           plyr_1.8.1           proto_0.3-10
[34] RColorBrewer_1.0-5   reshape2_1.4         rpart_4.1-8
[37] RSQLite_1.0.0        scales_0.2.4         sendmailR_1.2-1
[40] splines_3.1.0        stringr_0.6.2        survival_2.37-7
[43] tools_3.1.0          XML_3.98-1.1         xtable_1.7-4
[46] XVector_0.6.0

deseq2 paired • 11k views
4
Entering edit mode
@mikelove
Last seen 1 hour ago
United States

hi Hanni,

There is a trick from the edgeR user guide (under the section "Comparisons Both Between and Within Subjects") for building a R formula for such an experiment. For the above column data, add a new column for culture:

dds$culture.nested = factor(rep(rep(1:5,each=2),2)) Make sure that unstimulated is the base level of treatment: dds$treatment <- relevel(dds\$treatment, "unstimulated")

so now the column data looks like:

     tissue culture    treatment culture.nested
1  tissueX       1 unstimulated              1
2  tissueX       1   stimulated              1
3  tissueX       2 unstimulated              2
4  tissueX       2   stimulated              2
5  tissueX       3 unstimulated              3
6  tissueX       3   stimulated              3
7  tissueX       4 unstimulated              4
8  tissueX       4   stimulated              4
9  tissueX       5 unstimulated              5
10 tissueX       5   stimulated              5
11 tissueY       6 unstimulated              1
12 tissueY       6   stimulated              1
13 tissueY       7 unstimulated              2
14 tissueY       7   stimulated              2
15 tissueY       8 unstimulated              3
16 tissueY       8   stimulated              3
17 tissueY       9 unstimulated              4
18 tissueY       9   stimulated              4
19 tissueY      10 unstimulated              5
20 tissueY      10   stimulated              5

Then use the design:

design(dds) = ~ tissue + tissue:culture.nested + tissue:treatment

To make the extraction of results simpler, specify a standard model matrix:

dds = DESeq(dds, modelMatrixType="standard")

To test for genes which respond to treatment effect in tissue X:

results(dds, name="tissuetissueX.treatmentstimulated")

To test for genes which respond to treatment effect in tissue Y:

results(dds, name="tissuetissueY.treatmentstimulated")

To test for genes where the treatment effect differs between the two tissues:

results(dds, contrast=list("tissuetissueY.treatmentstimulated", "tissuetissueX.treatmentstimulated"))
0
Entering edit mode

Hi Michael,

Although this is a 5 years old post, it is still a relevant one. I have been reading all your answers regarding paired multifactor design. I hope, I understood most of it. What I am struggling to understand is if our goal is to test for genes which respond to treatment effect in different tissues then (tissue:treatment) interaction term is the one to go for, why will we need (tissue:culture.nested) interaction term?

Nurun

0
Entering edit mode

The difference between those terms is subtle and I’d recommend working with a statistician that can walk you through it.

0
Entering edit mode
hwillenbrock ▴ 10
@hwillenbrock-6937
Last seen 6.8 years ago
Denmark

Hi Michael,

Thanks!

That works, but won't this make the model assume that culture 1 and 6 are the same, as if the tissue samples where taken from the same person (which they are not)?

0
Entering edit mode

It's fine because this difference is absorbed in the main tissue effect. For the comparisons you want to make, the interactions are all you need.

0
Entering edit mode

Hi,

I believe I have an analogous issue, however I am still not sure what the best way to proceed is and would be grateful for your help.

In my multi-factor design  variables are:
treatment (before/after) / patient / disease (A/B) / response (Y/N)

For each patient I have one sample taken before treatment and one after (therefore with the same value for response and for disease).

1) The first question is to compare expression after/before treatment
I have been told by collaborators that the two diseases are expected to behave very differently and it does not make much sense to look at effects considering both together.
Thus, to keep things simple I created two different datasets and used design=~patient+ treatment in each of the them in order to extract the global effect of treatment across all patients, at the same time controlling for patient specific effects.
However, doing so I am ignoring the 'response' variable.

2) The second question is comparing responders (response=Y) vs non-responders, separately before and after treatment
One Naive solution would be to concentrate on only a subportion of the data (e.g. before treatment) and use design=~response
Instead, using all data with a design in which both patient and response variables are present gives problems with matrix rank.

Could you advise on the best design to answer both 1) and 2) ? And further clarify how the 'trick' to avoid problems with matrix rank should be implemented?

Thanks a lot!

0
Entering edit mode

1) If you want to compare expression before and after treatment, your current approach seems fine.

2) You can use the trick above to see if responders have a different treatment effect, controlling for patient differences. I would continue to split by disease to make things easier. Then follow the steps above, but for your analysis, "culture" is "patient" and "tissue" is "response".

Your design should look like ~ response + response:patient.nested + response:treatment. If you have different number of patients for response = Y and N then you'll have to do some extra work to create your own model matrix using model.matrix() and then to remove columns of the model matrix which have all zeros. See this other thread for how to do that:

C: DESeq2: Extracting interaction term

Run DESeq(dds, betaPrior=FALSE) or if you needed to create your own model matrix because the numbers of responders and non responders were different, run DESeq(dds, full=mm, betaPrior=FALSE). You will have the treatment effect in nonresponders and the treatment effect responders as the two interaction terms of treatment and response. These can be extracted with 'name' like above and contrasted with contrast=list(...) like above. The contrast tests for which genes the treatment effect was different for responders vs non-responders.

0
Entering edit mode

Hi Michael,

I would be extremely grateful if you could verify whether I implemented your suggestion correctly and clarify a few doubts regarding the extraction of results

This is how my projectinfo table looks like

sample treatment patient response patient.nested
2   control       1        Y             1
7   treated       1        Y             1
6   control       2        Y             2
10   treated      2        Y             2
5   control       3        Y             3
11   treated      3        Y             3
1   control       4        N             1
8   treated       4        N             1
4   control       5        N             2
9   treated       5        N             2
3   control       6        N             3
12   treated      6        N             3

Differential expression analysis:

> data <-DESeqDataSetFromMatrix(countData = countsTable,
+                                 colData = projectinfo,
+                                 design = ~response+response:patient.nested+response:treatment)
> diffexpr<-DESeq(data, test="Wald", betaPrior=FALSE)
> resultsNames(diffexpr)
[1] "Intercept"                  "response_Y_vs_N"            "responseN.patient.nested2"
[4] "responseY.patient.nested2"  "responseN.patient.nested3"  "responseY.patient.nested3"
[7] "responseN.treatmenttreated" "responseY.treatmenttreated"

The lists I want to extract

1a) genes with differential expression before/after treatment in responders
i.e. resp=Y&treatment=treated vs resp=Y&treatment=control
Is it correct that this is given by:
results(diffexpr,  name="responseY.treatmenttreated")  ?

1b) and analogously for non-responders (resp=N&treatment=treated vs resp=N&treatment=control) would be

results(diffexpr,  name="responseN.treatmenttreated")

2) genes with differential expression responders/non-responders separately for treated and controls, i.e.

2a) resp=Y&treatment=treated vs resp=N&treatment=treated
and
2b) resp=Y&treatment=control vs resp=N&treatment=control

It is not clear to me whether I can get these out of

results(diffexpr, contrast=list("responseY.treatmenttreated", "responseN.treatmenttreated"))

I realize that these are questions related to linear modeling and not specific to your tool, so I need to improve my knowledge on the topic and hope to find a nearby statistician who can help.

Thanks a lot again for your help!

Fulvia

0
Entering edit mode

1a) and 1b) yes

2a) and 2b) you cannot compare directly across response like this with a model which has fixed effects for controlling patient differences. How can you know if what you see if a change for response=Y and not specific to patient 1,2,3 in response? This you might want to discuss with a statistician if it's hard to see why you can't make this comparison.

The comparison I was referring to, where you contrast the two interaction terms, answers the question: is the treatment effect different for response=Y vs response=N.

0
Entering edit mode

Thanks for all the feedback! I have in the meantime been thinking about this and the importance of controlling for patient-to-patient differences.

0
Entering edit mode
@alicebdennis-7066
Last seen 5.3 years ago
Switzerland

Hi all,

I have a similar problem/nested design, but I think my problem is how to read in colData in the first place!

Is it possible to post your successful code in full here?

Thank you,

Alice

0
Entering edit mode

hi Alice,

You can post new questions with the "Ask question" button at the top. Your question here is posted as an answer to Hanni's question.

See this section of our DESeq2 workflow (the workflow goes at a slower pace than the vignette, explaining some base R functions as well):

"Typically, we have a table with experimental metadata for our samples. For your own project, you might create such a comma-separated value (CSV) file using a text editor or spreadsheet software such as Excel.

We load this file with read.csv. The parentheses around the last line are used to print the result in addition to storing it to the sampleTable object."

http://bioconductor.org/help/workflows/rnaseqGene/

If you didn't use summarizeOverlaps or htseq-count (i.e. you've already constructed a count matrix) you'll have to make sure that the column names of the counts matrix line up with the row names of the imported sample table.

0
Entering edit mode
hwillenbrock ▴ 10
@hwillenbrock-6937
Last seen 6.8 years ago
Denmark

Hi Michael,

One of my donor samples might be an outlier and I tried removing the samples from the analysis. This means that I now have the following samples:

    tissue culture    treatment culture.nested
1  tissueX       1 unstimulated              1
2  tissueX       1   stimulated              1
3  tissueX       2 unstimulated              2
4  tissueX       2   stimulated              2
5  tissueX       3 unstimulated              3
6  tissueX       3   stimulated              3
7  tissueX       4 unstimulated              4
8  tissueX       4   stimulated              4
9  tissueX       5 unstimulated              5
10 tissueX       5   stimulated              5
11 tissueY       6 unstimulated              1
12 tissueY       6   stimulated              1
13 tissueY       7 unstimulated              2
14 tissueY       7   stimulated              2
15 tissueY       8 unstimulated              3
16 tissueY       8   stimulated              3
19 tissueY      10 unstimulated              5
20 tissueY      10   stimulated              5

Now, when I try to do this:

design(dds) = ~ tissue + tissue:culture.nested + tissue:treatment

I get the error: "the model matrix is not full rank, so the model cannot be fit as specified."

Is there any way to do the analysis, still exploiting the paired design without having to remove the one of the tissueX cultures as well?

Thanks,

Hanni

0
Entering edit mode

hi Hanni,

Unfortunately, the base R function tries to form an interaction between tissue and all the non-base levels of culture.nested, leaving you with an interaction effect tissueY.culture.nested4 and no samples to fit that effect. I am planning to build more support for custom design matrices this devel cycle, for the cases like these where handing a formula to model.matrix() is not sufficient. You can use the workaround I propose at the link below, where you edit the design matrix to remove the column which cannot be fit:

C: DESeq2 paired and multi factor comparison

format: plain


text:
# Isotropic crystal and energy band

1. ### hokhani

300
Since all the directions are equivalent in an isotropic crystal, can we deduce that the energy band is exactly spherical?

2. ### DrDu

4,161
The problem is that there are no isotropic crystals. Even cubic crystals are only symmetric with respect to four and threefold rotations about some special axes. However this is sufficient to render second order (but not higher order) tensors isotropic.

3. ### hokhani

300
Do you mean that using the isotropic approximation we in fact disregard the crystal?

4. ### DrDu

4,161
Yes, indeed.

1 person likes this.


text:
# Just like 1980

There have been (conservatively) 16 “Moore generations” since 1980 – that is to say computing speed should have increased by $2^{16}$ – approximately 4 million times.

But computation requirements grow to fill the computing power available and BINSIC now runs at about the same speed, when executing the same code, as my ZX80 did when hacking its way through BASIC programs 32 years ago.

It’s a strange sight to see, because it is so evocative of that period.

The code runs slowly largely, I think, because of the way that I have chosen to handle `GOTO` statements (not that I could think of any other way) – essentially the whole or at least a substantial part, of the program is reparsed every time a `GOTO` is issued.

Anyway, it adds to the sense of realism!


text:
How to deal with the relative phase from the evaluation of the not interacting qubit?

Every qubit is a realization of some two-level quantum system with $$\left| 0 \right\rangle$$ and $$\left| 1 \right\rangle$$ states. These states have their energies $$E_0$$ and $$E_1$$ respectively. By solving the Schrödinger equation for the two-level quantum system one can show time evaluation of the $$\left| 0 \right\rangle$$, $$\left| 1 \right\rangle$$ and $$\left| + \right\rangle$$ states:

$$\left| \psi_0(t) \right\rangle = e^{-i \frac{E_0}{\hbar} t} \left| 0 \right\rangle \\ \left| \psi_1(t) \right\rangle = e^{-i \frac{E_1}{\hbar} t} \left| 1 \right\rangle \\ \left| \psi_+(t) \right\rangle = \frac{1}{\sqrt{2}}\left( e^{-i \frac{E_0}{\hbar} t} \left| 0 \right\rangle + e^{-i \frac{E_1}{\hbar} t} \left| 1 \right\rangle \right)= \frac{e^{-i \frac{E_0}{\hbar} t}}{\sqrt{2}} \left(\left| 0 \right\rangle + e^{-i \frac{E_1 - E_0}{\hbar} t} \left| 1 \right\rangle \right)$$

My question is about the relative phase $$\varphi = e^{-i \frac{E_1 - E_0}{\hbar} t}$$. This $$\phi$$ acts as a phase gate and changes the quantum state. Is this relative phase important and how in real quantum computers we should deal with it? One problem that I see is the following circuit:

After the measurement, the state should be $$\left| 0 \right\rangle$$, because $$HH = I$$. But If $$t$$ is such that the $$\varphi = \pi$$ then we will measure $$\left| 1 \right\rangle$$. We can forget about it if $$\varphi(t) << \pi$$ in the whole computational time. Otherwise, we should make sure that the time between gates should be $$2\pi$$.

Is this a problem and if yes how IBM, Google or etc. deal with this problem when constructing the circuits? How I understand the frequency $$(E_1 - E_0)/\hbar$$ value is important and for example, the ibmq_armonk qubit's frequency is equal to ~4.97428 GHz (qiskit_textbook), so should we take into account this value when we play with ibmq_armonk by using OpenPulse?

The could be a problem, but it depends on how you're realising your qubits. Some realisations are configured so that $$E_0=E_1$$, and then there's no problem.

There is (at least from the theoretical perspective) a simple fix: if you're supposed to be waiting a time $$t$$, then, instead:

• wait time $$t/2$$
• apply bit flip
• wait time $$t/2$$
• apply bit flip.

This exchanges the roles of 0 and 1 for half the time so that they end up with equal phases.

Alternatively, you just take this into account and apply compensating phase gates inside your circuit.

• Thanks for the answer. Can you mention an example of a qubit realization that has $E_0 = E_1$? – Davit Khachatryan Feb 28 '20 at 11:38
• IIRC, something like the hyperfine ground states of a Rubidium atom (e.g. trapped in an optical lattice) are a very good approximation to that because the timescales of the gates are much shorter than the timescale due to $1/(E_0-E_1)$. – DaftWullie Feb 28 '20 at 11:58
• You may wish to save your answer tick in case someone can answer your specific question about the ibm hardware, which is not something I kno anything about. – DaftWullie Feb 28 '20 at 11:59
• Thanks a lot. I will keep the answer tick for a while, but still, your answer was very helpful. – Davit Khachatryan Feb 28 '20 at 12:11


text:
# Not for the weak-hearted

Discrete Mathematics Level 4

Find the coefficient of $$x^{11}$$ in the expansion of $$(1+x^2)^4(1+x^3)^7(1+x^4)^{12}$$.

×


text:
# Seismic anisotropy

Series Investigations in Geophysics Öz Yilmaz http://dx.doi.org/10.1190/1.9781560801580 ISBN 978-1-56080-094-1 SEG Online Store

A medium is anisotropic if its intrinsic elastic properties, measured at the same location, change with direction, and it is isotropic if its properties do not change with direction [1]. Most of seismic data analysis is based on the assumption that the subsurface behaves seismically isotropic. A case of anisotropy that is of interest in exploration seismology is the change in velocity with direction.

A medium is transversely isotropic if its elastic properties do not change in any direction perpendicular to an axis of symmetry. The usual meaning of seismic anisotropy is variation of seismic velocity, which itself depends on the elastic properties of the medium, with the direction in which it is measured [2]. Anisotropic variation of seismic velocity must not be confused with the source-receiver azimuthal variation of moveout velocity for a dipping reflector in an isotropic medium (equation 7-3).

There are two cases of seismic anisotropy that we shall review in this section, both of which are special cases of transverse isotropy. For convenience, consider a horizontally layered earth model. First, is the vertical transverse isotropy, otherwise simply referred to as transverse isotropy, for which velocities do not vary from one lateral direction to another, but vary from one direction to another on a vertical plane that coincides with a given lateral direction. Horizontal bedding and fracturing parallel to the bedding produce transverse isotropy. Second, is the horizontal transverse isotropy, otherwise known as azimuthal anisotropy, for which velocities vary from one lateral direction to another. Fracturing in a direction other than the bedding direction gives rise to azimuthal anisotropy.

Transverse isotropy stems from the fact that, as a result of a depositional process, velocities in a layer are different in the direction of bedding and the direction perpendicular to the bedding. Specifically, fine layering of isotropic beds that constitute the depositional unit give rise to an anisotropy with its axis of symmetry perpendicular to the bedding plane. Azimuthal anisotropy can stem from the fact that, as a result of a tectonic process, rock material associated with a layer may have different rigidity in different azimuthal directions. A specific example is secondary fracturing in rocks whereby the velocity in the fracture direction is higher than the velocity in the orthogonal direction, again, giving rise to an anisotropy with its axis of symmetry parallel to the bedding plane.

Refer to Figure 11.7-1 to review the physical aspects of wave propagation in an anisotropic medium. The directional change of the velocity is illustrated by the skewed ellipse in Figure 11.7-1a, with the fast velocity in the direction of the major axis and the slow velocity in the direction of the minor axis. Now consider Huygens’ secondary sources situated along the wavefront A, say at time t. This wavefront actually coincides with an exploding reflector dipping at an angle ϕ. Because of the anisotropic behavior of the propagation medium, the Huygens’ sources do not emanate semicircular wavefronts; instead, the wavefronts are skewed in the direction of the fast velocity. These skewed wavefronts will form the plane wavefront B at a later time t + Δt. While the energy was transmitted along the raypath SP at group velocity, the wavefront that represents a constant phase actually traveled from position A to B along TP normal to the wavefront at phase velocity. Because the group velocity is associated with the raypath, it is sometimes referred to as the ray velocity. Similarly, because the phase velocity is associated with the wavefront, it is sometimes referred to as the wavefront velocity. Note that the wavefront angle θ associated with the phase velocity is different from the ray angle ϕ associated with the group velocity. Only if the medium were isotropic, Huygens’ secondary sources would produce semicircular wavefronts and the phase angle θ would coincide with the ray angle ϕ.

Note from Figure 11.7-1 that the zero-offset raypath SP does not make a right angle with the reflector; thus, in the case of anisotropy, the zero-offset ray is not a normal-incident ray as would be in the case of isotropy. This behavior can be better explained by sketching the isotropic and anisotropic wavefronts as shown in Figure 11.7-1b. Specifically, in the case of isotropic medium, the wavefront emanating from a point P is circular and the zero-offset ray is normal-incident to the reflector. Whereas, in the case of an anisotropic medium, the wavefront is skewed and the zero-offset ray impinges on the reflector at a non-normal incidence angle.

An example of how velocity changes with direction in an anisotropic medium is shown in Figure 11.7-2 [3]. Based on velocity characteristics of the Green River shale [4], this polar plot of group velocities shows that, as for most rock types, the P-wave velocity is nonelliptical, whereas the SH-wave velocity behaves elliptically anisotropic. Also note that the horizontal P-wave velocity is greater than the vertical P-wave velocity.

We shall associate seismic anisotropy primarily with seismic velocities. Aside from anisotropic velocity analysis, this means that we will need to review those processes that are intimately influenced by velocity anisotropy, such as dip-moveout correction, migration and AVO analysis.

The generalized form of Hooke’s law, which is the foundation of linear elastic theory, states that each stress component can be expressed as a linear combination of all the strain components [5]. Hooke’s law is based on the assumption that elastic deformations in solids are infinitesimally small. The stiffness matrix, otherwise known as the elastic modulus matrix [4], that relates the stress components to the strain components is (Section L.1)

 ${\displaystyle \{c_{ij}\}={\begin{pmatrix}c_{11}&c_{12}&c_{13}&c_{14}&c_{15}&c_{16}\\c_{21}&c_{22}&c_{23}&c_{24}&c_{25}&c_{26}\\c_{31}&c_{32}&c_{33}&c_{34}&c_{35}&c_{36}\\c_{41}&c_{42}&c_{43}&c_{44}&c_{45}&c_{46}\\c_{51}&c_{52}&c_{53}&c_{54}&c_{55}&c_{56}\\c_{61}&c_{62}&c_{63}&c_{64}&c_{65}&c_{66}\\\end{pmatrix}}.}$ (82a)

The elements of the stiffness matrix are the elastic constants of an elastic solid. Since this matrix is symmetric, cij = cji, there are 21 independent constants for an elastic medium.

For an isotropic solid, the elastic behavior of which is independent of direction at a point within the solid, the number of independent elastic constants is only two, known as Lamé’s constants, λ and μ. Consequently, the stiffness matrix given by equation (82a) reduces to the special form

 ${\displaystyle \{c_{ij}\}={\begin{pmatrix}\lambda +2\mu &\lambda &\lambda &0&0&0\\\lambda &\lambda +2\mu &\lambda &0&0&0\\\lambda &\lambda &\lambda +2\mu &0&0&0\\0&0&0&2\mu &0&0\\0&0&0&0&2\mu &0\\0&0&0&0&0&2\mu \\\end{pmatrix}}.}$ (82b)

To describe the P- and S-wave propagation in isotropic solids, only two elastic parameters are needed — Lamé’s constants, λ and μ. Other elastic parameters — Young’s modulus E, Poisson’s ratio σ and bulk modulus κ, and the P- and S-wave velocities can all be expressed in terms of λ and μ (Figure 11.0-1 and Appendix L.1).

For a transversely isotropic solid, the elastic behavior of which is the same in two orthogonal directions but different in the third direction, the number of independent constants is five [5] [4] [2]. For the more specific case of vertically transverse isotropy (VTI), which has a vertical symmetry axis, the five independent elastic constants are c11, c13, c33, c44, and c66 [4].

To explicitly describe the effect of anisotropy in wave propagation, Thomsen [4] has elegantly redefined the five elastic constants for the VTI media — the vertical P- and S-wave velocities, α0 and β0, in the vertical direction,

 ${\displaystyle \alpha _{0}={\sqrt {\frac {c_{33}}{\rho }}},}$ (83a)

 ${\displaystyle \beta _{0}={\sqrt {\frac {c_{44}}{\rho }}},}$ (83b)

and three constants that describe the degree of anisotropy, ε, γ, and δ, in terms of the five constants c11, c13, c33, c44, and c66

 ${\displaystyle \varepsilon ={\frac {c_{11}-c_{33}}{2c_{33}}},}$ (83c)

 ${\displaystyle \delta ={\frac {(c_{13}+c_{44})^{2}-(c_{33}-c_{44})^{2}}{2c_{33}(c_{33}-c_{44})}},}$ (83d)

and

 ${\displaystyle \gamma ={\frac {c_{66}-c_{44}}{2c_{44}}}.}$ (83e)

For most sedimentary rocks, the parameters ε, δ, and γ are of the same order of magnitude and usually much less than 0.2 [4]. In fact, the Thomsen parameters given by equations (83c,83d,83e) relate to the case of weak anisotropy described by small values (≪ 1) of ε, γ, and δ. While applications of anisotropy in seismic data analysis are primarily based on the assumption of weak anisotropy, these parameters still are useful to describe the general case of transverse isotropy.

The special case of δ = ε is known as elliptical anisotropy [7]. The ellipticity is associated with the shape of the wavefront expanding from a point source. Albeit its underlying theory is simpler than the general theory of anisotropy, elliptical anisotropy occurs in nature only rarely. Figure 11.7-3 shows a crossplot of the two Thomsen parameters, ε and δ, for various types of sedimentary and crystalline rocks based on field and laboratory studies [4]. Note that a few of the rock samples closely satisfy the ellipticity condition, ε = δ. Also, for most rock types, the Thomsen parameters, ε and δ, are positive constants less than 0.2; thus, the supporting evidence for weak anisotropy theory.

## References

1. Winterstein, 1990, Winterstein, D. F., 1990, Velocity anisotropy terminology for geophysicists: Geophysics, 55, 1070–1088.
2. Sheriff, 1991, Sheriff, R. E., 1991, Encyclopedic dictionary of exploration geophysics: Soc. Expl. Geophys.
3. Uren et al., 1990a, Uren, N. F., Gardner, G. H. F. and McDonald, J. A., 1990a, The migrator’s equation for anisotropic media: Geophysics, 55, 1429–1434.
4. Thomsen, 1986, Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, 1954–1966.
5. Officer, 1958, Officer, C. B., 1958, Introduction to the theory of sound transmission with application to the ocean: McGraw-Hill Book Co.
6. Uren et al., 1990b, Uren, N. F., Gardner, G. H. F. and McDonald, J. A., 1990b, Dip moveout in anisotropic media: Geophysics, 55, 863–867.
7. Daley and Hron, 1979, Daley, P. F. and Hron, F., 1979, Reflection and transmission coefficients for seismic waves in ellipsoidally anisotropic media: Geophysics, 44, 27–38.


text:
# Computing Extreme Eigenvalues of Large Scale Hankel Tensors

@article{Chen2016ComputingEE,
title={Computing Extreme Eigenvalues of Large Scale Hankel Tensors},
author={Yannan Chen and Liqun Qi and Qun Wang},
journal={Journal of Scientific Computing},
year={2016},
volume={68},
pages={716-738}
}
• Published 28 April 2015
• Mathematics
• Journal of Scientific Computing
Large scale tensors, including large scale Hankel tensors, have many applications in science and engineering. In this paper, we propose an inexact curvilinear search optimization method to compute Z- and H-eigenvalues of mth order n dimensional Hankel tensors, where n is large. Owing to the fast Fourier transform, the computational cost of each iteration of the new method is about $$\mathcal {O}(mn\log (mn))$$O(mnlog(mn)). Using the Cayley transform, we obtain an effective curvilinear search…
Computing Eigenvalues of Large Scale Sparse Tensors Arising from a Hypergraph
• Computer Science
SIAM J. Sci. Comput.
• 2016
A first-order optimization algorithm for computing H- and Z-eigenvalues of these large scale sparse tensors (CEST) and it is proved that the sequence of iterates generated by CEST converges to an eigenvector of the tensor.
Computing Z-Eigenvalue Of Hankel Tensors
• Computer Science
2019 16th International Computer Conference on Wavelet Active Media Technology and Information Processing
• 2019
This paper proposes a conjugate gradient algorithm with simple line search (CGLS) for solving Z-eigenvalue of Hankel tensors and analyzes the convergence of the CGLS algorithm and adjusts suitable parameters in the algorithm.
On computing minimal H-eigenvalue of sign-structured tensors
• Computer Science, Mathematics
• 2017
A sum-of-squares (SOS) algorithm for computing the minimal H-eigenvalues of tensors with some sign structures called extended essentially non negative tensors (EEN-tensors), which includes nonnegative tensors as a subclass.
Computing the generalized eigenvalues of weakly symmetric tensors
• Computer Science, Mathematics
Comput. Optim. Appl.
• 2017
Two convergent gradient projection methods for computing some generalized eigenvalues of weakly symmetric tensors and a spectral gradient projection method which is born from the combination of the BB method with the gradient projections method is superior to the GEAP, AG and AGP methods.
Tensor convolutions and Hankel tensors
• Mathematics, Geology
• 2017
Let A be an mth order n-dimensional tensor, where m, n are some positive integers and N:= m(n−1). Then A is called a Hankel tensor associated with a vector v ∈ ℝN+1 if Aσ = vk for each k = 0, 1,...,N
Finding the maximum eigenvalue of a class of tensors with applications in copositivity test and hypergraphs
• Mathematics, Computer Science
• 2015
A structured semi-definite program algorithm for computing the maximum eigenvalue of a class of tensors with sign structure called W-tensors, based on a new structured sums-of-squares decomposition result for a nonnegative homogeneous polynomial induced by a \$W-tensor.
A semidefinite program approach for computing the maximum eigenvalue of a class of structured tensors and its applications in hypergraphs and copositivity test
• Computer Science
Numer. Linear Algebra Appl.
• 2018
This paper introduces a new class of structured tensors called W‐tensors, which not only extends the well‐studied nonnegative tensors by allowing negative entries but also covers several important tensors arising naturally from spectral hypergraph theory.
The Decomposition Algorithm of Hankel Tensors
• Computer Science
2018 IEEE International Conference on Signal Processing, Communications and Computing (ICSPCC)
• 2018
An algorithm is discussed that can compute the decomposition expression and the rank of the Hankel tensor by iterative method and the efficiency and accuracy of proposed method can be demonstrated through numerical examples.
• Computer Science
Comput. Optim. Appl.
• 2016
In this paper, an adaptive gradient (AG) method is presented for generalized tensor eigenpairs and global convergence and linear convergence rate are established under some suitable conditions.
A trust region algorithm for computing extreme eigenvalues of tensors
• Computer Science
Numerical Algebra, Control & Optimization
• 2020
A trust region algorithm for the spherically constrained optimization of tensor eigenvalue problem to a nonlinear optimization with a spherical constraint is proposed in this paper.

## References

SHOWING 1-10 OF 64 REFERENCES
Fast Hankel tensor–vector product and its application to exponential data fitting
• Mathematics
Numer. Linear Algebra Appl.
• 2015
A fast algorithm for Hankel tensor–vector products is obtained by embedding a Hankel Tensor into a larger anti-circulant tensor, and the computational complexity is about O(m2nlogmn) for a square Hanko tensor of order m and dimension n.
Shifted Power Method for Computing Tensor Eigenpairs
• Computer Science
SIAM J. Matrix Anal. Appl.
• 2011
A shifted symmetric higher-order power method (SS-HOPM), which it is shown is guaranteed to converge to a tensor eigenpair, and a fixed point analysis is used to characterize exactly which eigenpairs can and cannot be found by the method.
A sequential subspace projection method for extreme Z-eigenvalues of supersymmetric tensors
• Mathematics
Numer. Linear Algebra Appl.
• 2015
The main idea of SSPM is to form a 2-dimensional subspace at the current point and then solve the original optimization problem in the subspace and the globalization strategy of random phase can be easily incorporated into S SPM, which promotes the ability to find extreme Z-eigenvalues.
Computing Tensor Eigenvalues via Homotopy Methods
• Computer Science, Mathematics
SIAM J. Matrix Anal. Appl.
• 2016
An upper bound is derived for the number of equivalence classes of generalized tensor eigenpairs using mixed volume and an algorithm is introduced that combines a heuristic approach and a Newton homotopy method to extract real generalized eigpairs from the found complex generalized eigen pairs.
Efficient MATLAB Computations with Sparse and Factored Tensors
• Computer Science
SIAM J. Sci. Comput.
• 2007
This paper considers how specially structured tensors allow for efficient storage and computation, and proposes storing sparse tensors using coordinate format and describes the computational efficiency of this scheme for various mathematical operations, including those typical to tensor decomposition algorithms.
On the Best Rank-1 Approximation of Higher-Order Supersymmetric Tensors
• Mathematics
SIAM J. Matrix Anal. Appl.
• 2002
It is shown that a symmetric version of the above method converges under assumptions of convexity (or concavity) for the functional induced by the tensor in question, assumptions that are very often satisfied in practical applications.
Z-eigenvalue methods for a global polynomial optimization problem
• Computer Science
Math. Program.
• 2009
This paper proposes some Z-eigenvalue methods for solving the problem of the best rank-one approximation to higher order tensors, and proposes a direct orthogonal transformation Z- eigenvalue method for this problem in the case of order three and dimension three.
Hankel Tensors: Associated Hankel Matrices and Vandermonde Decomposition
Hankel tensors arise from applications such as signal processing. In this paper, we make an initial study on Hankel tensors. For each Hankel tensor, we associate it with a Hankel matrix and a higher
Finding the extreme Z-eigenvalues of tensors via a sequential semidefinite programming method
• Mathematics
Numer. Linear Algebra Appl.
• 2013
An approximation method is proposed to solve the TCLP by using a sequence of semidefinite programming problems, and it gives a numerical algorithm to compute the extreme Z-eigenvalue of an even order tensor with dimension larger than three, which improves the literature.


text:
# How to apply fourier transform to $0.5^n u(n)$

I'm working in a signals class for continuous signals, and we have this problem shown above. I have tried using this function $f_1*f_2 = F_1 * F_2$, where I'm assuming this means multiplication of two functions is equal to the convolution of their fourier transforms. I'm using $f_1 = 0.5^n$ and $f_2 = u(n)$.

So I can calculate the fourier transorm of $u(n)$ fine. It is $\pi\delta(\omega) + 1/(j\omega))$. However, I cannot for the life of me figure out $0.5^n$. I tried to put it into the fourier transform integral integral of$(0.5^t)/(e^{j \omega t})dt$ from negative infinity to infinity, but I end up with $0.5t/e^{jw}$, and when evaluated from negative infinity to infinity, I end up with $\infty$ as my answer, unless of course the integration is wrong.

Therefore, either the answer is $\infty * \pi\delta(w) + 1/(j\omega)$, which when convoluted would equal just the second function..? OR am I going about this problem completely wrong?

I think you'd better not split the functions. If you want to calculate the fourier transform of $0.5^t u(t)$, then you might put all that into the integral. Therefore, your integral limits will become $0$ to $+\infty$ because of $u(t)$. And your integral result should not be infinity because your function $0.5^t$ goes to zero.

Solving the integral:

$\large\int_0^\infty e^{(-jwt)} . 0.5^t dt$

the inside part could be rewrite:

$\huge\frac{e^{(-jwt)}}{2^t} => (\frac{e^{(-jw)}}{2})^t$

You could name $\large\frac{e^{(-jw)}}{2} = u$ so you get integral of $\large u^t dt = \frac{u^t}{ln(u)} + C$

Then you get (replacing u):

$\huge\frac{(\frac{e^{-jw}}{2})^t}{ ln(e^{-jw}/2)}$

When t goes to $+\infty$ you get zero. When it goes to zero you get:

$\large\frac{1}{ln(e^{-jw}/2)}$

which left us:

$\large\frac{1}{-jw - ln(2)}$

I'm not sure I did not get lost in calculations but I believe its correct.

• You're completely right. I got mixed up after separating the functions. Thank you! – nathpilland Nov 25 '13 at 2:12
• You are welcome. Good luck – Felipe_Ribas Nov 25 '13 at 3:26


text:
### Comments:

Very nice paper! The main idea: if the distance between two embedded nodes has to be large then one should not waste energy in order to make this distance a bit smaller. Trying to make large distances a bit smaller is actually what Laplacian embedding is doing with its quadratic objective. One should rather focus on small distances and give up on the pairs of embedded nodes that are far away. ### Embedding complex networks: This method could be applied to social networks (or complex networks in general): social networks are known to have "weak ties" (roughly speaking, weak ties are links between people that do not have many friends in common) that may perturb laplacian embedding. ### Asymptotic complexity and scalability in practice: The suggested algorithm to solve the optimization does not seem very scalable. The asymptotic running time complexity is not given. ### Concerns about the optimization in dimension higher than 1. The example in 2-D. In equations (9-11), what if we change the constraints into: $||x||^2+||y||^2=1$, $e^Tx=0$ and $e^Ty=0$? In n-D, denoting $r_i$ the vector associated to node $i$, that would give $\sum_ir_i^2=1$ and $\sum_ir_i=0$. This seems to be a more relevant optimization to me: minimize the distance between nodes in 2-D without setting constraints for each dimension, but rather global constraints. If the objective is the one of the laplacian embedding, is it true that, in the optimal solution, $y$ will be proportional to $x$? ($x$ and $y$ will both be proportional to the eigenvector associated with the smallest eigenvalue of the Laplacian?) If so, is it still true if the objective is another one? Like the Cauchy embedding objective for instance? ### Others: - How to tune the parameter $\sigma$?
> - How to tune the parameter $σ$? I don't see where $σ$ is defined, did I miss something?
Note that $\sum_i r_i = 0$ means that vector $r$ is orthogonal to vector $v_0 = (1,1,\ldots\,1,1,1)$. The latter corresponds to the smallest eigen-value $\lambda_0 = 0$ of [Laplacian matrix](https://en.wikipedia.org/wiki/Laplacian_matrix). Adding constraint $r \perp v_0$ to the optimization problem drives us to seek the second smallest eigenpair. You may also read very related and well-written [lecture](https://ocw.mit.edu/courses/mathematics/18-409-topics-in-theoretical-computer-science-an-algorithmists-toolkit-fall-2009/lecture-notes/MIT18_409F09_scribe3.pdf) about Rayleigh quotients, min-max theorem (Courant–Fischer–Wey). Gives a good background in the field of eigen-based solutions for optimization problems. In 2016, during the preparation of [my talk](http://www.complexnetworks.fr/temporal-density-of-complex-networks-and-ego-community-dynamics/) titled "Temporal density of complex networks and ego-community dynamics" ([slides](https://kirgizov.link/talks/temporal-density-4-july-2016.pdf), french). I have read stuff about graph optimization problems. Below I present my personal (perhaps a little bit random selection) of related papers and books * [Four Cheeger-type Inequalities for Graph Partitioning Algorithms](https://www.math.ucsd.edu/~fan/wp/heaticcm.pdf), Fan Chung, 2007 * [Normalized cuts and image segmentation](https://people.eecs.berkeley.edu/~malik/papers/SM-ncut.pdf), Jianbo Shi and Jitendra Malik, 2000 * [Eigenvalues of graphs](https://web.cs.elte.hu/~lovasz/eigenvals-x.pdf), László Lovász, 2007 * [Spectres de graphes](https://www-fourier.ujf-grenoble.fr/~ycolver/All-Articles/98a.pdf) de Yves Colin de Verdière, 1998
Read the paper, add your comments…

### Comments:

Read the paper, add your comments…

### Comments:

## Great overview of graph embedding methods ### The encoder, decoder and loss function formalism is very well designed! ### Similarity matrix: In Table 1, maybe there are methods dedicated to the case where the similarity matrix is required to be a positive semidefinite (SDP) matrix. Indeed, (i) there are plenty of such relevant SDP matrices (such as "$D+A$", "$A*A$" or "the cosine similarity" which are all SDP) and (ii) the optimal embedding minimizing the L2 loss function is easy to find using the eigenvectors associated with the largest eigenvalues. It thus makes this kind of method very handy. This kind of [kernel PCA](https://en.wikipedia.org/wiki/Kernel_principal_component_analysis)-like method does not work if the similarity matrix (somehow the kernel here) is not SDP. ### Shallow embedding VS autoencoder-based embedding: Limitations of "shallow embedding methods" are highlighted on page 9 and it seems that autoencoder-based methods answer some of them. But this part is not totally clear to me. From what I understand, the major difference between the two kinds of methods is that: - for the methods referred to as "shallow embedding methods", if you want to reconstruct the entry $s(i,j)$ of the input similarity matrix, then you need the learned feature vectors of both nodes $i$ and $j$ and then the reconstruction of $s(i,j)$ is often simple, e.g. the scalar product between the two feature vectors (cf. equation (5)) or sometimes a softmax (cf. equation (10)) is added on top of the scalar product; - while in the case of autoencoder-based embedding methods, you only need the learned feature vector of a single node $i$ to reconstruct the full column $s_i$ of the input similarity matrix and the reconstruction is more elaborated. We can notice that no parameters are shared between nodes in the decoder for the shallow methods, while some parameters can be shared in the case of autoencoder-based methods. There thus can be fewer parameters in autoencoder-based methods and thus it could be in principle faster and have some regularisation. In addition, if a previously unseen node arrives, an autoencoder-based method could in principle still give a feature vector for that node, while it does not seem to be the case for shallow methods (they are referred to as "inherently transductive" methods). "Shallow embedding also fails to leverage node attributes during encoding". I think that it is possible to adapt both shallow methods and autoencoder-based methods to take into account this additional information on node attributes. I didn't get this point. As noted in the survey, it seems that current autoencoder-based methods (SDNE and DNGR) are not scalable and are inherently transductive: "the input dimension to the autoencoder is fixed at $|V|$, which can be extremely costly and even intractable for graphs with millions of nodes. In addition, the structure and size of the autoencoder is fixed, so SDNE and DNGR are strictly transductive and cannot cope with evolving graphs, nor can they generalize across graphs." The scalability issue has also been noted in [VERSE](https://papers-gamma.link/paper/48/): "Works such as [12, 48] investigate deep learning approaches for graph embeddings. Their results amount to complex models that require elaborate parameter tuning and computationally expensive optimization, leading to time and space complexities unsuitable for large graph" where [12] and [48] are SDNE and DNGR. ### Neighborhood aggregation and convolutional encoders: This part seems important. In particular, Algorithm 1 seems important as it is referred to in subsequent parts. However, I found this section very hard to understand. ### Incorporating task-specific supervision: Even though this section is short, I find it very interesting. ### Typos: - "they found the the more complex aggregators" - In equation (24) $h_j^{k-1}$ - End of page 20. What is $\Epsilon$ in $O(|\Epsilon|)$? - "(e.g. , using..." - "[8] J. Bruna, W. Zaremba, and Y. Szlam, A.and LeCun. Spectral networks and locally connected networks on graphs..."
Read the paper, add your comments…


text:
# What is the jackknife?

According to Wikipedia, the jackknife technique “was developed by Maurice Quenouille from 1949 and refined in 1956″, and was so named by John Tukey in 1958 (who also expanded on the technique). Quenouille’s 1956 paper “Notes on Bias in Estimation” gives a clear explanation of the intuition behind the jackknife estimate.

In the general statistical set-up, we are interested in estimating some parameter $\theta$, using observations $x_1, \dots, x_n$. Let $\hat\theta_n$ denote our estimate for $\theta$ based on these $n$ observations. (For concreteness, one example you could consider is $\theta = \sigma^2$, i.e. variance of the distribution, and $\hat\theta_n = \frac{1}{n}\sum_{i=1}^n (x_i - \bar{x})^2$.)

There are several different properties we might want our estimate $\hat\theta_n$ to have, e.g. efficiency, sufficiency, consistency, unbiasedness. The maximum likelihood estimate (MLE) has the first 3 properties but not the 4th (except in rare situations). If we have an MLE (or some other consistent estimate), how can we adjust it to obtain an estimate that is less biased?

Under some conditions, we can perform a Taylor expansion to get

\begin{aligned} \mathbb{E}[\hat\theta_n - \theta] = \frac{a_1}{n} + \frac{a_2}{n^2} + \frac{a_3}{n^3} + \dots \end{aligned}

If this is the case, then if we define $\theta' = n\hat\theta_n - (n-1)\hat\theta_{n-1}$, the equation above gives

\begin{aligned} \mathbb{E}[\theta'] &= n \left( \theta + \frac{a_1}{n} + \frac{a_2}{n^2} + \frac{a_3}{n^3} + \dots \right) - (n-1)\left( \theta + \frac{a_1}{n-1} + \frac{a_2}{(n-1)^2} + \frac{a_3}{(n-1)^3} + \dots \right) \\ &= \left( n\theta + a_1 + \frac{a_2}{n} + \frac{a_3}{n^2} + \dots \right) - \left( (n-1)\theta + a_1 + \frac{a_2}{n-1} + \frac{a_3}{(n-1)^2} + \dots \right) \\ &= \theta + \frac{a_2}{n(n-1)} + \dots \\ &= \theta + O\left( \frac{1}{n^2} \right), \end{aligned}

so $\theta'$ is biased to order $O\left( \frac{1}{n^2} \right)$, while our original estimate $\hat\theta_n$ was biased to order $O\left( \frac{1}{n} \right)$.

In defining $\theta'$, we used $\hat\theta_{n-1}$ which was based on $x_1, \dots, x_{n-1}$. There is no a priori reason why we should this particular set of $n-1$ points: we could have easily used $\{ x_2, \dots, x_n \}$, $\{ x_1, x_3, x_4, \dots, x_n \}$, or any other set of $n-1$ points instead. In fact, Quenouille argues that in order to minimize the loss in efficiency, one should use all possible sets of $n-1$ observations.

If we define $\hat\theta_{(-i)}$ as the estimate of $\theta$ based on all the observations except $x_i$, then replacing $\hat\theta_{n-1}$ by the mean of these $\hat\theta_{(-i)}$‘s (which we denote by $\tilde\theta$) yields the jackknife estimate:

\begin{aligned} \hat\theta_{jack} = n\hat\theta_n - (n-1) \tilde\theta = n\hat\theta_n - (n-1) \cdot \frac{1}{n}\sum_{i=1}^n \hat\theta_{(-i)}. \end{aligned}

Instead of bias estimation, Tukey saw the utility of the jackknife for robust interval estimation. Arguing that the $\hat\theta_{(-i)}$‘s could be treated as approximately i.i.d. random variables, the statistic

$\dfrac{\sqrt{n}(\tilde\theta - \theta)}{\sqrt{\frac{1}{n-1} \sum_{i=1}^n (\hat\theta_{(-i)} - \tilde\theta)^2} }$

should have an approximate $t$ distribution with $n-1$ degrees of freedom. Thus, we can make use of this pivotal statistic to do interval estimation.

References:

1. Miller, R. (1974). The Jackknife–A Review.
2. Quenouille, M. H. (1956). Notes on Bias in Estimation.
3. Wikipedia. Jackknife resampling.

# Beamer: inserting section slides before each section

When creating a presentation in LaTeX using Beamer, the \tableofcontents automatically creates a table of contents for you based on the \section, \subsection and \subsection tags. For example, in one of my documents,

\begin{frame}{Outline of talk}
\tableofcontents
\end{frame}

creates this slide:

Sometimes, in addition to a table of contents slide at the beginning of the presentation, I want a slide at the beginning of each section to remind the audience of where we are in the presentation. We can do this by adding the following lines of code before the \begin{document} command:

\AtBeginSection[]
{
\begin{frame}
\tableofcontents[currentsection]
\end{frame}
}

This will create slides just like \tableofcontents, except that all but the current section will be grayed out. For example, the second section will start with the following slide:

If you want this behavior for subsections as well, you can add the following code:

\AtBeginSubsection[]
{
\begin{frame}
\tableofcontents[currentsection,currentsubsection]
\end{frame}
}

References:

# The math behind quantile regression

In this previous post, we introduced quantile regression and demonstrated how to run it in R. In this post, we dig into the math behind quantile regression.

Actually, before we get into quantile regression, let’s talk about quantiles. Given some observations $y_1, y_2, \dots, y_n \in \mathbb{R}$, how do we define the 50th quantile, better known as the median? In high school, we are taught that if we line up all our observations in ascending order, the median is the element right in the middle if $n$ is odd, or the sum of the two middle elements if $n$ is even. Quantiles can be defined in a similar manner: for $\tau \in (0,1)$, the $\tau$th quantile is a value such that at least $\tau$ fraction of the observations are less than or equal to it.

This is all well and good, but it’s a bit hard to generalize this notion easily. Instead, we can define the $\tau$th quantile as the solution to a minimization problem; the “algorithm” for finding the $\tau$th quantile in the previous paragraph can be viewed as a way to solve this minimization problem.

Here is the minimization problem: if we define $\rho_\tau(y) = y (\tau - 1\{ y < 0\})$, then for a set of observations $y_1, \dots, y_n \in \mathbb{R}$, we defined its $\tau$th quantile as

\begin{aligned} \hat{q}_\tau &= \underset{q \in \mathbb{R}}{\text{argmin}} \quad \sum_{i=1}^n \rho_\tau (y_i - q) \\ &= \underset{q \in \mathbb{R}}{\text{argmin}} \quad \sum_{i=1}^n (y_i - q)(\tau - 1\{ y_i < q \}) \\ &= \underset{q \in \mathbb{R}}{\text{argmin}} \quad \left[ (\tau - 1) \sum_{y_i < q} (y_i - q) + \tau \sum_{y_i \geq q} (y_i - q) \right]. \quad (1) \end{aligned}

When $\tau = 1/2$, we see that the above reduces to

\begin{aligned} \hat{q}_{0.5} &= \underset{q \in \mathbb{R}}{\text{argmin}} \quad \left[ -\frac{1}{2} \sum_{y_i < q} (y_i - q) + \frac{1}{2} \sum_{y_i \geq q} (y_i - q) \right] \\ &= \underset{q \in \mathbb{R}}{\text{argmin}} \quad \sum_{i=1}^n |y_i - q|. \end{aligned}

Here is some hand-wavy reasoning for why the sample $\tau$th quantile minimizes the expression in $(1)$. Assume that a proportion $p$ of the terms are in the first summation $\displaystyle\sum_{y_i < q} (y_i - q)$, and that $1-p$ of the terms are in the second summation $\displaystyle\sum_{y_i \geq q} (y_i - q)$. If we increase $q$ by some little $\epsilon$ so that the number of terms in each summation doesn’t change, then overall the expression increases by

\begin{aligned} np(1-\tau)\epsilon - n(1-p)\tau \epsilon = n\epsilon (p - \tau). \end{aligned}

If $p < \tau$, then we can decrease the objective function value in $(1)$ by making $q$ larger. If $p > \tau$, then we can decrease the objective function value by making $q$ smaller. The only place where we can’t have any improvement is when $p = \tau$, i.e. $q$ chosen such that $\tau$ of the observations are less than it.

Let’s extend this minimization problem idea to quantile regression. For each response value $y_i$, we have a corresponding estimate $\hat{y}_i$ We can interpret the setting above as trying to minimize

$\displaystyle\sum_{i=1}^n \rho_\tau (y_i - \hat{y}_i),$

where the $\hat{y}_i$‘s are all forced to be equal to one value, $q$. In quantile regression, we are forced instead to have $\hat{y}_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}$, where $x_{i1}, \dots, x_{ip}$ are the values for the $i$th observation for the $p$ features on hand. Thus, we can rewrite the minimization problem as

$\underset{\beta_0, \dots, \beta_p}{\text{minimize}} \quad \displaystyle\sum_{i=1}^n \rho_\tau (y_i - \beta_0 - \beta_1 x_{i1} - \dots - \beta_p x_{ip}),$

or in vector notation (and incorporating the intercept as a feature),

$\underset{\beta}{\text{minimize}} \quad \displaystyle\sum_{i=1}^n \rho_\tau (y_i - x_i^T \beta).$

Notice that this is very similar to ordinary least squares regression. In that case, instead of applying $\rho_\tau$ to each $y_i - x_i^T \beta$, we apply $x \mapsto \frac{1}{2}x^2$.

References:

1. Wikipedia. Quantile regression.

# Quantile regression in R

Quantile regression: what is it?

Let $y \in \mathbb{R}$ be some response variable of interest, and let $x \in \mathbb{R}^p$ be a vector of features or predictors that we want to use to model the response. In linear regression, we are trying to estimate the conditional mean function, $\mathbb{E}[y \mid x]$, by a linear combination of the features.

While the conditional mean function is often what we want to model, sometimes we may want to model something else. On a recent episode of the Linear Digressions podcast, Katie and Ben talked about a situation in Uber where it might make sense to model a conditional quantile function.

Let’s make this more concrete. Say Uber came up with a new algorithm for dispatching drivers and we are interested in how this algorithm fares in terms of wait times for consumers. A simple (linear regression) model for this is

$\mathbb{E}[wait\_time \mid algorithm] = a + b \cdot algorithm,$

where $algorithm = 1$ if the new algorithm was used to dispatch a driver, and $algorithm = 0$ if the previous algorithm was used. From this model, we can say that under the old algorithm, mean wait time was $a$, but under the new algorithm, mean wait time is $a + b$. So if $b < 0$, I would infer that my new algorithm is "doing better".

But is it really? What if the new algorithm improves wait times for 90% of customers by 1 min, but makes the wait times for the remaining 10% longer by 5 min? Overall, I would see a decrease in mean wait time, but things got significantly worse for a segment of my population. What if that 10% whose wait times became 5 minutes longer were already having the longest wait times to begin with? That seems like a bad situation to have, but our earlier model would not pick it up.

One way to pick up such situations is to model conditional quantile functions instead. That is, trying to estimate the mean of $y$ given the features $x$, let’s trying to estimate a quantile of $y$ given the features $x$. In our example above, instead of trying to estimate the mean wait time, we could estimate the 95th quantile wait time to catch anything going wrong out in the tails of the distribution.

Another example where estimating conditional quantiles is useful is in growth charts. All of you have probably seen one of these charts below in a doctor’s office before. Each line in the growth chart represents some quantile for length/weight given the person’s age. We track our children’s length and weight on this chart to see if they are growing normally or not.

WHO growth chart for boys. The lines represent conditional quantile functions for different quantiles: of length given age and weight given age.

Quantile regression is a regression method for estimating these conditional quantile functions. Just as linear regression estimates the conditional mean function as a linear combination of the predictors, quantile regression estimates the conditional quantile function as a linear combination of the predictors.

Quantile regression in R

We can perform quantile regression in R easily with the quantreg package. I will demonstrate how to use it on the mtcars dataset. (For more details on the quantreg package, you can read the package’s vignette here.)

Let’s load our packages and data:

library(quantreg)
data(mtcars)

We can perform quantile regression using the rq function. We can specify a tau option which tells rq which conditional quantile we want. The default value for tau is 0.5 which corresponds to median regression. Below, we fit a quantile regression of miles per gallon vs. car weight:

rqfit <- rq(mpg ~ wt, data = mtcars)
rqfit
# Call:
#     rq(formula = mpg ~ wt, data = mtcars)
#
# Coefficients:
#     (Intercept)          wt
# 34.232237   -4.539474
#
# Degrees of freedom: 32 total; 30 residual

Printing the fitted object to the console gives some rudimentary information on the regression fit. We can use the summary function to get more information (just like we do for the lm function):

summary(rqfit)
# Call: rq(formula = mpg ~ wt, data = mtcars)
#
# tau: [1] 0.5
#
# Coefficients:
#     coefficients lower bd upper bd
# (Intercept) 34.23224     32.25029 39.74085
# wt          -4.53947     -6.47553 -4.16390

In the table above, the lower bd and upper bd columns represent the endpoints of confidence intervals for the model coefficients. There are a number of ways for these confidence intervals to be computed; this can be specified using the seoption when invoking the summary function. The default value is se="rank", with the other options being “iid”, “nid”, “ker”, “boot” and “BLB” (type ?summary.rq for details).

This next block of code plots the quantile regression line in blue and the linear regression line in red:

plot(mpg ~ wt, data = mtcars, pch = 16, main = "mpg ~ wt")
abline(lm(mpg ~ wt, data = mtcars), col = "red", lty = 2)
abline(rq(mpg ~ wt, data = mtcars), col = "blue", lty = 2)
legend("topright", legend = c("lm", "rq"), col = c("red", "blue"), lty = 2)

Median regression (i.e. 50th quantile regression) is sometimes preferred to linear regression because it is “robust to outliers”. The next plot illustrates this. We add two outliers to the data (colored in orange) and see how it affects our regressions. The dotted lines are the fits for the original data, while the solid lines are for the data with outliers. As before, red is for linear regression while blue is for quantile regression. See how the linear regression fit shifts a fair amount compared to the median regression fit (which barely moves!)?

y <- c(mtcars$mpg, 40, 36) x <- c(mtcars$wt, 5, 4)
plot(y ~ x, pch = 16, main = "mpg ~ wt")
points(c(5, 4), c(40, 36), pch = 16, col = "dark orange")
abline(lm(mpg ~ wt, data = mtcars), col = "red", lty = 2)
abline(lm(y ~ x), col = "red")
abline(rq(mpg ~ wt, data = mtcars), col = "blue", lty = 2)
abline(rq(y ~ x), col = "blue")

As I mentioned before, the tau option tells rq which conditional quantile we want. What is interesting is that we can set tau to be a vector and rq will give us the fits for all those quantiles:

multi_rqfit <- rq(mpg ~ wt, data = mtcars, tau = seq(0, 1, by = 0.1))
multi_rqfit
# Call:
#     rq(formula = mpg ~ wt, tau = seq(0, 1, by = 0.1), data = mtcars)
#
# Coefficients:
#     tau= 0.0  tau= 0.1 tau= 0.2  tau= 0.3  tau= 0.4  tau= 0.5  tau= 0.6  tau= 0.7  tau= 0.8  tau= 0.9
# (Intercept)  27.6875 37.509794  37.2768 34.876712 34.891176 34.232237 36.734405 38.497404 38.879916 44.391047
# wt           -3.7500 -6.494845  -6.2400 -5.205479 -4.852941 -4.539474 -5.016077 -5.351887 -5.250722 -6.266786
# tau= 1.0
# (Intercept) 44.781558
# wt          -5.627981
#
# Degrees of freedom: 32 total; 30 residual

I don’t think there’s a clean one-liner to plot all the quantile fits at once. The following loop is one possible solution:

# plotting different quantiles
colors <- c("#ffe6e6", "#ffcccc", "#ff9999", "#ff6666", "#ff3333",
"#ff0000", "#cc0000", "#b30000", "#800000", "#4d0000", "#000000")
plot(mpg ~ wt, data = mtcars, pch = 16, main = "mpg ~ wt")
for (j in 1:ncol(multi_rqfit\$coefficients)) {
abline(coef(multi_rqfit)[, j], col = colors[j])
}

rq exhibits interesting behavior if tau is set to some value outside of $[0,1]$: it gives solutions for ALL values of tau in $[0,1]$! (It turns out that there are a finite number of them.)

One final note: rq also supports more complicated formula.

# no intercept
rq(mpg ~ wt - 1, data = mtcars)
# Call:
#     rq(formula = mpg ~ wt - 1, data = mtcars)
#
# Coefficients:
#     wt
# 4.993498
#
# Degrees of freedom: 32 total; 31 residual

rq(mpg ~ wt + cyl, data = mtcars)
# Call:
#     rq(formula = mpg ~ wt + cyl, data = mtcars)
#
# Coefficients:
#     (Intercept)          wt         cyl
# 38.871429   -2.678571   -1.742857
#
# Degrees of freedom: 32 total; 29 residual

# model with interactions
rq(mpg ~ wt * cyl, data = mtcars)
# Call:
#     rq(formula = mpg ~ wt * cyl, data = mtcars)
#
# Coefficients:
#     (Intercept)          wt         cyl      wt:cyl
# 51.6650791  -7.2496674  -3.4759342   0.5960682
#
# Degrees of freedom: 32 total; 28 residual

# fit on all other predictors
rq(mpg ~ ., data = mtcars)
# Call:
#     rq(formula = mpg ~ ., data = mtcars)
#
# Coefficients:
#     (Intercept)         cyl        disp          hp        drat          wt        qsec          vs          am
# 7.61900909  0.64658010  0.02330302 -0.03149324  0.78058219 -4.56231735  0.57610198  1.58875367  1.33175736
# gear        carb
# 2.37453015 -0.33136465
#
# Degrees of freedom: 32 total; 21 residual

References:

1. University of Virginia Library. Getting started with quantile regression.
2. Wikipedia. Quantile regression.

# Simulated annealing for optimization

Simulated annealing is a probabilistic technique for approximating the global optimum of a given objective function. Because it does not guarantee a global optimum, it is known as a metaheuristic.

According to Wikipedia, annealing is “a heat treatment that alters the physical and sometimes chemical properties of a material to increase its ductility and reduce its hardness, making it more workable. It involves heating a material above its recrystallization temperature, maintaining a suitable temperature for a suitable amount of time, and then cooling” (emphasis mine). As we will see, simulated annealing mimics this physical process in trying to find a function’s global optimum.

There are several versions of simulated annealing. Below, we give a high-level outline of a basic version for minimizing a function $f: \mathcal{X} \mapsto \mathbb{R}$. (Simulated annealing can be used for constrained optimization as well; for simplicity we limit our exposition to unconstrained optimization.) Let $i$ index the iteration number. Simulated annealing has a positive temperature parameter, denoted by $T_i$ at iteration $i$, whose evolution strongly influences the result.

Start at some initial point $x_0 \in \mathcal{X}$. For each iteration $i = 0, \dots$, if the stopping criterion is not satisfied:

1. Choose some random proposal point $y$ from a probability distribution. The lower the temperature, the more concentrated the probability distribution’s mass is around $x_i$.
2. With some acceptance probability, set $x_{i+1} = y$. If not, set $x_{i+1} = x_i$. The acceptance probability is a function of the function values $f(x_i)$ and $f(y)$, and the current temperature $T_i$.
3. Update the temperature parameter.

That’s all there is to it! The devil, of course, is in the details. Here are some examples of what each of the building blocks could be.

Stopping criterion:

• Terminate after a fixed number of iterations, or when the computational budget has been fully utilized.
• Terminate when the objective function does not improve much, e.g. $|f(x_i) - f(x_{i+1})| < 10^{-7}$.

Choosing a random proposal point: This depends very much on what the underlying state space looks like. A benefit of simulated annealing is that it can be applied to very general state spaces $\mathcal{X}$, even discrete ones, as long as there is a notion of distance between two points in the space. (In particular, first-order information about the space (i.e. gradients) is not needed.)

• If $\mathcal{X} = \mathbb{R}^d$, we could set $y = x_i + T \epsilon$, where $\epsilon \sim \mathcal{N}(0,I_d)$, or $\epsilon \sim U[-1, 1]^d$.
• The dependence on $T$ could be different as well, e.g. $y = x_i + \sqrt{T} \epsilon$.
• If $\mathcal{X}$ is the space of permutations on $\{ 1, 2, \dots, n\}$, then $y$ could be selected uniformly at random from the set of permutations which are at most $T$ transpositions from $x_i$.

In the proposals above, the smaller $T$ is, the closer $y$ is likely to be than $x_i$.

Acceptance probability: The probability of accepting $y$ as the next value in the sequence (i.e. as $x_{i+1}$) is a function which depends on $f(x_i)$, $f(x_{i+1})$ and (crucially) the temperature $T$. Typically if $f(y) < f(x_i)$, we would accept $y$ with probability 1, but this does not have to be the case. In general, the acceptance probability should decrease as $T$ decreases and as $f(y) - f(x)$ becomes more positive.

• This is the most popular acceptance function in practice: if $f(y) - f(x) < 0$, accept. If not, accept with probability $\dfrac{1}{1 + \exp \left( \frac{f(y) - f(x)}{T} \right)}$.

Updating the temperature parameter: We just require $T_i \rightarrow 0$ as $i \rightarrow \infty$.

• $T_i = T_0 \cdot (0.95)^i$.
• $T_i = T_0 / i$.
• $T_i = T_0 / \log i$.

References:

1. Henderson et al. The theory and practice of simulated annealing.
2. Mathworks. How simulated annealing works.
3. Mathworks. Simulated annealing terminology.
4. Wikipedia. Simulated annealing.


text:
### Economic costs of context collapse

One motivation for my recent studies on information flow economics is that I’m interested in what the economic costs are when information flows across the boundaries of specific markets.

For example, there is a folk theory of why it’s important to have data protection laws in certain domains. Health care, for example. The idea is that it’s essential to have health care providers maintain the confidentiality of their patients because if they didn’t then (a) the patients could face harm due to this information getting into the wrong hands, such as those considering them for employment, and (b) this would disincentivize patients from seeking treatment, which causes them other harms.

In general, a good approximation of general expectations of data privacy is that data should not be used for purposes besides those for which the data subjects have consented. Something like this was encoded in the 1973 Fair Information Practices, for example. A more modern take on this from contextual integrity (Nissenbaum, 2004) argues that privacy is maintained when information flows appropriately with respect to the purposes of its context.

A widely acknowledged phenomenon in social media, context collapse (Marwick and boyd, 2011; Davis and Jurgenson, 2014), is when multiple social contexts in which a person is involved begin to interfere with each other because members of those contexts use the same porous information medium. Awkwardness and sometimes worse can ensue. These are some of the major ways the world has become aware of what a problem the Internet is for privacy.

I’d like to propose that an economic version of context collapse happens when different markets interfere with each other through network-enabled information flow. The bogeyman of Big Brother through Big Data, the company or government that has managed to collect data about everything about you in order to infer everything else about you, has as much to do with the ways information is being used in cross-purposed ways as it has to do with the quantity or scope of data collection.

It would be nice to get a more formal grip on the problem. Since we’ve already used it as an example, let’s try to model the case where health information is disclosed (or not) to a potential employer. We already have the building blocks for this case in our model of expertise markets and our model of labor markets.

There are now two uncertain variables of interest. First, let’s consider a variety of health treatments $J$ such that $m = \vert J \vert$. The distribution of health conditions in society is distributed such that the utility of a random person $i$ receiving a treatment $j$ is $w_{i,j}$. Utility for one treatment is not independent from utility from another. So in general $\vec{w} \sim W$, meaning a person’s utility for all treatments is sampled from an underlying distribution $W$.

There is also the uncertain variable of how effective somebody will be at a job they are interested in. We’ll say this is distributed according to $X$, and that a person’s aptitude for the job is $x_i \sim X$.

We will also say that $W$ and $X$ are not independent from each other. In this model, there are certain health conditions that are disabling with respect to a job, and this has an effect on expected performance.

I must note here that I am not taking any position on whether or not employers should take disabilities into account when hiring people. I don’t even know for sure the consequences of this model yet. You could imagine this scenario taking place in a country which does not have the Americans with Disabilities Act and other legislation that affects situations like this.

As per the models that we are drawing from, let’s suppose that normal people don’t know how much they will benefit from different medical treatments; $i$ doesn’t know $\vec{w}_i$. They may or may not know $x_i$ (I don’t yet know if this matters). What $i$ does know is their symptoms, $y_i \sim Y$.

Let’s say person $x_i$ goes to the doctor, reporting $y_i$, on the expectation that the doctor will prescribe them treatment $\hat{j}$ that maximizes their welfare:

$\hat j = arg \max_{j \in J} E[X_j \vert y]$

Now comes the tricky part. Let’s say the doctor is corrupt and willing to sell the medical records of her patients to her patient’s potential employers. By assumption $y_i$ reveals information both about $w_i$ and $x_i$. We know from our earlier study that information about $x_i$ is indeed valuable to the employer. There must be some price (at least within our neoclassical framework) that the employer is willing to pay the corrupt doctor for information about patient symptoms.

We also know that having potential employers know more about your aptitudes is good for highly qualified applicants and bad for not as qualified applicants. The more information employers know about you, the more likely they will be able to tell if you are worth hiring.

The upshot is that there may be some patients who are more than happy to have their medical records sold off to their potential employers because those particular symptoms are correlated with high job performance. These will be attracted to systems that share their information across medical and employment purposes.

But for those with symptoms correlated with lower job performance, there is now a trickier decision. If doctors are corrupt, it may be that they choose not to reveal their symptoms accurately (or at all) because this information might hurt their chances of employment.

A few more wrinkles here. Suppose it’s true the fewer people will go to corrupt doctors because they suspect or know that information will leak to their employers. If there are people who suspect or know that the information that leaks to their employers will reflect on them favorably, that creates a selection effect on who goes to the doctor. This means that the information that $i$ has gone to the doctor, or not, is a signal employers can use to discriminate between potential applicants. So to some extent the harms of the corrupt doctors fall on the less able even if they opt out of health care. They can’t opt out entirely of the secondary information effects.

We can also add the possibility that not all doctors are corrupt. Only some are. But if it’s unknown which doctors are corrupt, the possibility of corruption still affects the strategies of patients/employees in a similar way, only now in expectation. Just as in the Akerlof market for lemons, a few corrupt doctors ruins the market.

I have not made these arguments mathematically specific. I leave that to a later date. But for now I’d like to draw some tentative conclusions about what mandating the protection of health information, as in HIPAA, means for the welfare outcomes in this model.

If doctors are prohibited from selling information to employers, then the two markets do not interfere with each other. Doctors can solicit symptoms in a way that optimizes benefits to all patients. Employers can make informed choices about potential candidates through an independent process. The latter will serve to select more promising applicants from less promising applicants.

But if doctors can sell health information to employers, several things change.

• Employers will benefit from information about employee health and offer to pay doctors for the information.
• Some doctors will discretely do so.
• The possibility of corrupt doctors will scare off those patients who are afraid their symptoms will reveal a lack of job aptitude.
• These patients no longer receive treatment.
• This reduces the demand for doctors, shrinking the health care market.
• The most able will continue to see doctors. If their information is shared with employers, they will be more likely to be hired.
• Employers may take having medical records available to be bought from corrupt doctors as a signal that the patient is hiding something that would reveal poor aptitude.

In sum, without data protection laws, there are fewer people receiving beneficial treatment and fewer jobs for doctors providing beneficial treatment. Employers are able to make more advantageous decisions, and the most able employees are able to signal their aptitude through the corrupt health care system. Less able employees may wind up being identified anyway through their non-participation in the medical system. If that’s the case, they may wind up returning to doctors for treatment anyway, though they would need to have a way of paying for it besides employment.

That’s what this model says, anyway. The biggest surprise for me is the implication that data protection laws serve this interests of service providers by expanding their customer base. That is a point that is not made enough! Too often, the need for data protection laws is framed entirely in terms of the interests of the consumer. This is perhaps a politically weaker argument, because consumers are not united in their political interest (some consumers would be helped, not harmed, by weaker data protection).

References

Akerlof, G. A. (1970). The market for” lemons”: Quality uncertainty and the market mechanism. The quarterly journal of economics, 488-500.

Davis, J. L., & Jurgenson, N. (2014). Context collapse: theorizing context collusions and collisions. Information, Communication & Society, 17(4), 476-485.

Marwick, A. E., & Boyd, D. (2011). I tweet honestly, I tweet passionately: Twitter users, context collapse, and the imagined audience. New media & society, 13(1), 114-133.

Nissenbaum, H. (2004). Privacy as contextual integrity. Wash. L. Rev., 79, 119.