text:
## Wednesday, October 2, 2013

### First day of State of Food Insecurity (SOFI) 2013

The FAO flagship publication SOFI 2013 was release yesterday on the 1st of October, the publication is the most important report in monitoring the progress towards the 2015 Millenium Development Goal and ultimately eliminate hunger.

I was interest in how the people responded, so I scrapped some data from Twitter and previous work to carry out some basic analysis.

In total, there were 1,284 Tweets over the course of 24 hours. However, 86% were retweets which does not add much information for text mining and association analysis.

Nevertheless, lets have a look at what people are saying using some basic frequency and sentimental analysis. Clearly, hunger is the most mentioned words however there are words which are interesting and they have been colored according to theopinion lexicon of Hu and Liu. Negative feelings are red, positive feelings are blue while neutral words has remained grey.

The negative words like miss, bad and slow has pointed out that more efforts are required to achieve the MDG; furthermore people are suffer from chornic hunger. Nevertheless, the work of FAO has been recognised by good work and progress.

From the word cloud we can observe that it is impossible to use a single lexicon to generalize emotions. For example, fallen may usually used in a negative way but in this context it represents the positive movement of rate of undernourishment falling. A special lexicon may be required to examine the social feedback about SOFI.

Unfortunately I had problem obtaining the geocode from the API, otherwise it would be interesting to examine where these tweets originate.

Given that majority of the tweets are retweets of the official release message, little information can be obtained from these tweets. We hope a revision in a few weeks/months time will provide us with much more insight about the status of the publication.

The codes can be found at my Github repository .

## Tuesday, July 2, 2013

### Keys to win a basketball game

Offence sells tickets, defences wins games, rebounding wins championships - Pat Summitt

In fact, if we try to translate the game of basketball into mathematics (very naively) we have the following equation.

\begin{align} \text{Score advantage}_x &= \text{Possession}_x \times \text{Probability of Scoring}_x - \text{Possession}_y \times \text{Probability of Scoring}_y\\ &= [\text{out of bounds}_x + \text{rebound}_x] \times \text{Probability of Scoring}_x\\ &- [\text{out of bounds}_y + \text{rebound}_y] \times \text{Probability of Scoring}_y\\ &= [\text{out of bounds}_x + \text{rebound}_x] \times \text{Probability of Scoring}_x\\ &- [\text{out of bounds}_y + (\text{total rebound} - \text{rebound}_x)] \times \text{Probability of Scoring}_y \end{align}

From the equation we can see that a good scoring team has a high probability to score, while a good defensive team reduces the opponents ability and probability to score. Yet, we can see that the rebound appears on both end where a rebound for team x is one less rebound for team y. Which modifies the number of possession giving more chance to attack and a potentially larger score advantage!

Of course, rebounds can be offensive and defensive and they are a function of your and your opponents scoring ability. The events are inter-connected and should be modelled as a chained event. Maybe I'll do that when I recover from my hang-over.

## Sunday, March 24, 2013

### Tupper's self-referential formula

Can't remember where I first came across this equation but the Tupper's self referential equation, is a very interesting formula that when graphed in two dimension plane it reproduces the formula.

$\frac{1}{2} < \left\lfloor \bmod\left(\left\lfloor\frac{y}{17}\right\rfloor2^{-17\lfloor x\rfloor - \bmod(\lfloor y \rfloor, 17)}, 2\right)\right\rfloor$

I first thought this would be a quick 5 min exercise which turned into a 3 hour work, the obstacle was that the constant “k” used in the formula is an extremely big integer and can not be handled in R naturally.

After a little search Large integer in R and play around it seems that the gmp library seems to work well.

First you will need to install GMP (The GNU Multiple Precision Arithmetic Library) from http://gmplib.org/, then install the gmp library within R.

## Load the library after installing GMP
library(gmp)

## Define the constant k
k = as.bigz("960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719")

## The tupper's formula
tupper = function(x, y, k){
z1 = as.bigz(y + k)
z2 = as.bigq(z1/17)
z3 = 2^(-17 * x - as.bigz(z1%%17))
0.5 < floor(as.bigz(z2 * z3)%%2)
}

## The x and y axis
x = 0:105
y = 0:16

## Compute the matrix
a = matrix(0, nc = length(x), nr = length(y))
for(i in seq_along(x)){
a[ ,107 - i] = rev(tupper(x[i], y, k = k))
}

Here is the plot

image(t(a), col = c("white", "black"))

## Wednesday, March 20, 2013

### Violin plots and regional income distribution

While preparing my slides for statistical graphics, a plot really caught my eye when I was playing around with the data.

I started off by plotting the time seriesof GNI per capita by country, and as expected it got quite messy and incomprehensible.

## Download and manipulate the data
library(FAOSTAT)
raw.lst = getWDItoSYB(indicator = c("NY.GNP.PCAP.CD", "SP.POP.TOTL"))
raw.df = raw.lst[["entity"]]
traw.df = translateCountryCode(raw.df, from = "ISO2_WB_CODE", to = "UN_CODE")
mraw.df = merge(traw.df, FAOregionProfile[, c("UN_CODE", "UNSD_MACRO_REG")])
final.df = mraw.df[!is.na(mraw.df$UNSD_MACRO_REG), ] ## Simple ugly time series plot ggplot(data = final.df, aes(x = Year, y = NY.GNP.PCAP.CD)) + geom_line(aes(col = Country)) + labs(x = NULL, y = "GNI per capita")  So I decided to compute the weighted average by region to examine the regional trends. ## Compute regional aggregates based on UN M49 definition reg.df = aggRegion(aggVar = "NY.GNP.PCAP.CD", weightVar = "SP.POP.TOTL", data = traw.df, keepUnspecified = FALSE, aggMethod = "weighted.mean", relationDF = data.frame(UN_CODE = FAOregionProfile[, "UN_CODE"], REG_NAME = FAOregionProfile[, "UNSD_MACRO_REG"])) ## Plot regional aggregates ggplot(data = reg.df[!is.na(reg.df$NY.GNP.PCAP.CD), ],
aes(x = Year, y = NY.GNP.PCAP.CD)) +
geom_line(aes(col = REG_NAME)) +
labs(x = NULL, y = "GNI per capita", col = "")

I can now see the trend clearly, but there are two problems with this approach. First, the variability within region is vast and thus the weighted average or any summary statistic such as quantile can be misleading and it does not tell me what is going on within the regions. Secondly, since a minimum of 65% of the country must be present in order to compute the aggregation, no statistics was available prior to 1985.

While I was carrying out regional comparisons with box-plot and violin plot I thought why not plot them accross time as well! So here is the final graph:

## Time series violin plot
ggplot(data = final.df,
aes(x = as.character(Year), y = NY.GNP.PCAP.CD)) +
geom_violin() + scale_y_log10() +
facet_wrap(~UNSD_MACRO_REG, ncol = 1, scales = "free_y") +
scale_x_discrete(breaks = as.character((seq(1960, 2010, by = 10))),
labels = as.character((seq(1960, 2010, by = 10)))) +
labs(x = NULL, y = "GNI per capita")

Now I can compare the regions, but at the same time I can see the within region income distribution. It amazes me how the income distribution diverges in Europe and Oceania while America and Asia moves towards a bell shaped distribution. Growth in Africa appears to be slow, but there are several countries which are growing at a faster rate and pushing the tail of the distribution. Although some of the variability in the density may have resulted from independence of countries, nonetheless it is still infromative.

## Monday, February 11, 2013

### R package building automation

Title: R package building automation

Inspired by the post at http://giventhedata.blogspot.tw/2013/02/my-r-package-development-cheat-sheet.html. I have decided to publish my cheat script for package development as well.

Building package used to be a nightmare, filling in all those Rd files manually can cause some serious brain damage. Thanks to the roxygen2 team and knitr, building a package with nice documentation is merely running a script like what I have below. Just modify and change pkgName to the name of your package then your package should be ready in 5 minute.

## Build the package
## ---------------------------------------------------------------------

## Remove the folder if it exists
if(file.exists("./pkgName"))

## Build the package
package.skeleton("pkgName", code_files = paste("./Codes/",
dir("./Codes/", pattern = "\\.R$"), sep = ""), force = TRUE) ## Include the data dir.create("pkgName/data") file.copy(from = "./FAOcountryProfile.RData", to = "pkgName/data/", overwrite = TRUE) file.copy(from = "./FAOregionProfile.RData", to = "pkgName/data/", overwrite = TRUE) file.copy(from = "./FAOmetaTable.RData", to = "pkgName/data/", overwrite = TRUE) file.copy(from = "./DESCRIPTION", to = "pkgName/", overwrite = TRUE) ## Include Demo dir.create("pkgName/demo") file.copy(from = "./pkgNamedemo.R", to = "pkgName/demo/", overwrite = TRUE) cat("pkgNamedemo Demonstration for the pkgName package\n", file = "pkgName/demo/00Index") ## Use roxygen to build the documentation library(roxygen2) roxygenize("pkgName") ## Include vignette dir.create("./pkgName/inst/doc/") file.copy(from = "../Documentation/pkgName/pkgName.pdf", to = "./pkgName/inst/doc/", overwrite = TRUE) ## Create the vignette hack from (http://yihui.name/knitr/demo/vignette/) ## This is not required for R 3.0.0 cat("%\\VignetteIndexEntry{General Manual}\n\\documentclass{article}\n\\begin{document}\n\\end{document}", file = "./pkgName/inst/doc/pkgName.Rnw") ## Build and check the package system("R CMD INSTALL --build pkgName") system("R CMD build pkgName") system("Rcmd check pkgName") system("Rcmd check --as-cran pkgName")  PS: This is my first post using Rmd, the highlight for R is defintely much better than plain text. However, I think some twick is required for my blog theme. ## Tuesday, February 5, 2013 ### Relearn boxplot and label the outliers Despite the fact that box plot is used almost every where and taught at undergraduate statistic classes, I recently had to re-learn the box plot in order to know how to label the outliers. This stackoverflow post was where I found how the outliers and whiskers of the Tukey box plots are defined in R and ggplot2: In ggplot2, what do the end of the boxplot lines represent? and this post on how to label the outliers using base graphics. How to label all the outliers in a boxplot Since the use of ggplot2 is required for this task, I have written some basic hack code to label the outliers for ggplot2. Here are the codes:  ## Install the FAOSTAT package to obtain the data if(!is.element("FAOSTAT", .packages())) install.packages("FAOSTAT") library(FAOSTAT) ## Download data on Cassava production cp.lst = getFAOtoSYB(name = "cassava_production", domainCode = "QC", itemCode = 125, elementCode = 5510) ## Use the country level data, and take only data for 2011 and remove the NA's cp.df = cp.lst$entity[!is.na(cp.lst$entity$cassava_production) &
cp.lst$entity$Year == 2011, ]

## Merge with the country profile to obtain the country names for labelling
ccp.df = merge(cp.df, FAOcountryProfile[, c("FAOST_CODE", "ABBR_FAO_NAME")],
all.x = TRUE)

## Merge with the regional pofile to obtain the UNSD M49 macro region
## composition for multiple boxplot.
rcp.df = merge(ccp.df, FAOregionProfile[, c("FAOST_CODE", "UNSD_MACRO_REG")],
all.x = TRUE)

## Compute the quantile
qrcp.df = ddply(.data = rcp.df, .variables = .(UNSD_MACRO_REG), transform,
lQntl = quantile(cassava_production, probs = 0.25, na.rm = TRUE),
uQntl = quantile(cassava_production, probs = 0.75, na.rm = TRUE))

## Compute the lower and upper bound which defines the outlier
brcp.df = ddply(.data = qrcp.df, .variables = .(UNSD_MACRO_REG), transform,
lBound = lQntl - 1.5 * (uQntl - lQntl),
uBound = uQntl + 1.5 * (uQntl - lQntl))

## Remove the country names if it is within the bounds
with(brcp.df, {
brcp.df[cassava_production <= uBound &
cassava_production >= lBound, "ABBR_FAO_NAME"] <<- NA
})

## Plot the data
set.seed(587)
ggplot(data = brcp.df, aes(x = UNSD_MACRO_REG, y = cassava_production)) +
geom_boxplot(outlier.colour = NA) +
geom_text(aes(label = ABBR_FAO_NAME), size = 2,
position = position_jitter(width = 0.1)) +
labs(x = NULL, y = NULL, title = "Production of Cassava by region")

Here is the final product, to avoid over-plotting of texts I have used position_jitter. Which is not an elegant solution but I just can not find any algorithm that works well in general.

## Sunday, February 3, 2013

### A package for agricultural statistic: FAOSTAT

After 8 years of using R, today I finally become a contributor to the community and released my first package, FAOSTAT.

The package is designed to provide user with direct access to the FAOSTAT data base via R and to support the open data and methodology philosophy used in the Statistical Yearbook of the Food and Agricultural Organization.

The package can be downloaded and installed the same way as other package in a single line:
install.packages("FAOSTAT")

Or the latest development version from github:
 library(devtools)
install_github("FAOSTAT","mkao006")

A draft documentation can be found at:

vignette("FAOSTAT", package = "FAOSTAT")

I hope the package will save time for people working with agricultural statistics and more generally official statistics and to promote the use for research. All feedbacks and criticisms are welcome.

### Maize trade Part II: Comparison and analysis

Following my last post about the maize network, although interesting but is not very informative. What we are going to do today is to contrast the maize network with the wine trade network.

The choice why we have chose wine will become clear after the network and the analysis. Lets first have a look at the trade network of wine, again only the top 80% of the trade has been shown.

First, we can see that the wine network is much more inter-connected than the maize network. Secondly, the amount of reciprocal trading is much higher. Third, there does not appear to have a single trading of the size of the magnitude such as the export of maize from United state to Japan and Korea

The advantage of the network diagram not only reveals the connection, but the structure of the network, hence the nature of the market. It is clear that the market for maize and wine is different, but how, and why, is what we are interested.

Usage/consumption:
First of all, the usage of the two commodity are on the opposite end of the spectrum. Maize is a staple crop, necessary in many regional cuisine also used for livestock feed and bio-fuel. On the other hand, wine is consumed at a lesser extent and is considered as a luxury or enjoyment (although considered a necessity in Italy and several other countries).

Cost:
The usage of the commodity lead to the difference in cost and elasticity. Staples and biofuel need to be cheap to be viable and thus cost is an important consideration. In contrast, the preference and quality of wine plays a much bigger role than the actual price tag. This is where the trade agreement plays a significant role in promoting the trade of maize between the United States, Japan and Korea but a smaller element for wine trading.

Finally, the variety and the value added of the commodity characterizes the trade. Japan can import maize from any country in the world to use for consumption, feed or to produce bio-fuel without any considerable change or impact because the value add of maize is small; simply put, the product can be easily substituted. On the other hand, wine is dominated by several producers such as France, Italy and a handful of new world wine producers. Which is not easily substituted given the preferences. This is the contributing reason to the high number of trading partners, interconnected clusters and reciprocal trading; since a single trading partner can not satisfy the preferences of the consumer and even the French would like to have Italian wine occasionally.

This example illustrate how a simple network diagram can reveal a lot of interesting insight compare to analyzing just the quantity and value of trade, it also supports my view on how important exploratory analysis is.

stringsAsFactors = FALSE)

wineEx.df = subset(wine.df,subset = type == "Export",
c("reporting_country", "partner_country", "Wine"))

## Sort the data and take only the top 80% of the trade
wineEx.df = arrange(wineEx.df, desc(Wine))
wineEx.df$sWine = scale(wineEx.df$Wine, center = FALSE)
wineEx.df$cs = cumsum(wineEx.df$Wine)
wineEx.df[wineEx.df[, 1] == "China, Hong Kong SAR ", 1] = "China, Hong Kong SAR"
wineFinal.df = subset(wineEx.df, cs < tail(wineEx.df$cs, 1) * 0.8) ## Set edge and arrow size wineFinal.df$edgeSize = with(wineFinal.df, Wine/sum(Wine))
wineFinal.df$arrowSize = ifelse(wineFinal.df$edgeSize * 30 < 0.5 , 0.5,
wineFinal.df$edgeSize * 15) ## Create the network and plot wine.net = network(wineFinal.df[, 1:2]) plot(wine.net, displaylabels = TRUE, label.col = "steelblue", edge.lwd = c(wineFinal.df$edgeSize) * 100,


text:
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Are you an
Engineering professional?
Join Eng-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

# Lighting technologies?

## Lighting technologies?

(OP)
This is certainly not my field of expertise. So, I may be putting the question wrongly. Anyhow:

I have overheard several discussions about lighting lately. Mostly office and domestic lighting. There seems to be a lot of options these days. Incandescent is on its way down, LED on its way up, Field Effect Emission I don't know where it is going, low power/compact flourescent, induction(?) and probably also other technologies that I haven't come across. Mercury and sodium are for larger areas, I think. But if they are available also for office and homes, why not?

Anyone sitting on the whole truth? Or having a good link?

Gunnar Englund
www.gke.org
--------------------------------------
100 % recycled posting: Electrons, ideas, finger-tips have been used over and over again...

### RE: Lighting technologies?

Right now, I don't think there's that much in the way of choice for domestic interior lighting; it's pretty much CFL all the way.  For domestic exterior light, LEDs might be tolerable, for garden lighting, etc., although there's a lot of useless LED lights for that as well.  The MR11 class spot lighting with a minimum of 12 3528 SMD LEDs seem to work OK.  I still waiting on some wedge base lights for the other garden lights.

One major thing to note is that even the "warm" whites are still on the blue side of warm.

### RE: Lighting technologies?

Being an early adopter, I've repeatedly tried LED light bulbs. So far, they've always failed way too early. Based on that experience, I think we should give them a few more years to work the bugs out before trying again.

CFL bulbs are now the default. They're now available for about US$2 or less. Even the monster-sized 45w (150-200 watt equivalent) are as cheap as$10.

They seems to have sorted-out most of the early failures. And the retail prices are so low that we no longer wail-and-weep and hold funeral services for each failed CFL bulb.

### RE: Lighting technologies?

I have zero CFL bulbs installed and I intend to keep it that way. I have yet to see a decent color rendering or instant on CFL (I live in Michigan and it gets cold in some places). I do have LED bulbs on the other hand. Not all are that great as IRstuff pointed out but if you are willing to pay the extra price then the new Philips and GE bulbs are excellent (in my opinion). The Sylvania bulbs are OK but I don't care for the pinkish 60W equivalent. The 40W non-dimming equivalents were on clearance for $10 and I have those in a location that is difficult to get to. Harold SW2011 SP1.0 OPW2011 SP1.0 Win 7 Ultimate BOXX 8550 Xtreme Dual Xeon 5680 @4.2Ghz nVidia Quadro 4000 www.lumenflow.com ### RE: Lighting technologies? Skogs; I'm having a hard time fishing a question out of your post. The only one I see is: #### Quote: Mercury and sodium are for larger areas, I think. But if they are available also for office and homes, why not? Sodium and Mercury Vapor, while very efficient (hence their heavy application in area illumination), they have absolutely abysmal color rendering or CRI. Sodium is pretty much one dominant spectral line at 589.3nm bletch yellow. MV is much better but is still sports a horrid CRI. For years we've lit our living room with MH which has an excellent CRI over 90. That burned out a few weeks ago so I just built this: It's twenty four LEDs run by a current controller. That shot doesn't begin to convey the photon flux spewing out of that thing. It's about 4x the impact expressed here in this picture. Keith Cress kcress - http://www.flaminsystems.com ### RE: Lighting technologies? I've always noticed a significant noise level off the mercury or sodium lights (I'm not sure which). that is probably limiting for residential applications. John D ### RE: Lighting technologies? (OP) Smoked As I said, I am probably not putting the question the way it should be put. If there is any. The background is that a company has invited me to buy shares in a field emission type of lamp. I have told them that they probably won't be able to compete with LED and that their timing is bad. Their market window is closed and they should have moved a lot faster. They started the development back in 1996 when they took over a Russian patent. They have then consumed around 30 MUSD to develop that technology into a marketable device. They have, as of now, mot even been able to publish test data that are on par with LED fixtures. All they can say is that their DESIGN GOAL is better than what can be done with LEDs today. I think that company is a scam and I have told them so. They were not a scam from the beginning, but now when they probably understand that they can never reach their design goals, it has turned into a scam. So, I thought that it would be good to have your opinion. If you had mentioned Field Emission, that would perhaps have been a support for that company. But you didn't. And that supports my thinking. Gunnar Englund www.gke.org -------------------------------------- 100 % recycled posting: Electrons, ideas, finger-tips have been used over and over again... ### RE: Lighting technologies? Ah. There are a few oddballs out there and your assessment is correct. They are too little, tooooo late. It's my prediction none of them will amount to a hillah beans. There are so many aspects to lighting beyond just - light. Trust and perceptions of the buyers is really important. Going in to an architect with some off base lighting methodology is just not going to fly these days. LEDs have had a long slow acceptance and they're way more approachable than Field Emission or any of those other odd balls. I think you made a wise decision. Keith Cress kcress - http://www.flaminsystems.com ### RE: Lighting technologies? I wondered if that is who you were wondering about, Gunnar. FWIW - it got some mention somewhere here a few months back. I bought some shares of stock on spec., and a month or so later the share price bounced to roughly twice what I'd bought it for, and I sold it, making about$500 for my retirement plan.

It's since dropped back to its historical level as a penny stock ($.50 per share, US), so I bought another few$100 worth to hold, just as a "what the heck" move.  Represents way less than 1% of my retirement fund, so I can afford to lose it...and you can't win big if you don't bet on a few long shots...

The technology might find a niche, I figure, if they can price it somewhere close to LED's.  The directionality and color temperature should compare favorably to LED fixtures, which are too much like spotlights for my taste.  CFL's are going to become a nightmare soon, given that nobody (even here in tree-hugging Seattle) has really figured out how to keep them out of the landfills.  The total mercury content in them is low, but somebody soon is going to start complaining about it, and the real costs of CFL's (in terms of fuss and bother and media hype) is going to climb again.

Then again, I need to start looking into who is doing tallow rendering, as candles may start to compete again too.  Oh, and find out who owns those Japanese whaling trawlers, whale oil lamps used to be hot sellers...

### RE: Lighting technologies?

Not sure if the Field Effect Emision lamp is the same as the Electron Stimulated Luminescence lamp that I saw some info on last month from VU1.  Looked interesting.  This one was basically a small CRT without the bad sitcoms broadcast on it.

In the commercial US field of lighting, we're probably going to have linear and compact fluorescents for a long while for the majority of indoor general lighting.  Seeing LEDs used commercially for accent or decorative fixtures.

Residentially it's hard to get the masses away from thinking that all light sources have to be a medium screw base (yes I know bayonet is used in Europe).  When most people hear "fluorescent", they imagine the old "cool white" or "daylight" lamps of 30 years ago.  ("Daylight" was an awful color).

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Close Box

# Join Eng-Tips® Today!

Join your peers on the Internet's largest technical engineering professional community.
It's easy to join and it's free.

Here's Why Members Love Eng-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!

format: plain


text:
# Find 25th term of the arithmetic series 9, 15, 21, …

Find 25th term of the arithmetic series 9, 15, 21, …
You can still ask an expert for help

• Questions are typically answered in as fast as 30 minutes

Solve your problem for the price of one coffee

• Math expert for every subject
• Pay only if we can solve it

Kimberly Evans
Let $\left({x}_{n}\right)=\left(9,15,21,27,.......\right)$ be the arithmetic sequence.
The common difference is $d={x}_{n+1}-{x}_{n}=6$
The first term is a=9.
The general term is given by ${x}_{n}=a+\left(n-1\right)d$
$\therefore {x}_{25}=9+\left(24-1\right)\cdot 6$
$=147$


text:
# What is the pressure and density of a neutron star?

I am interested in finding pressure of neutron star! So. Please could any tell me how to choose central density for the inner and outer core of neutron Star. What numeric value should me in both core. Also what should be the polytropic index gamma for outer core.

• How can you have a "central density" for an "outer core"? Nov 8, 2016 at 7:37

The central density of a neutron star is unknown, but probably lies in the region of $\sim 10^{18}$ kg/m$^3$, depending on the details of the composition, the equation of state and of course the mass of the neutron star.
As an order of magnitude estimate, at these densities, the pressure approaches the maximum for a relativistically degenerate ideal fermion gas of $P \simeq \rho c^2/3 = 3\times10^{34}$ Pa.
By "outer core" I assume you mean the neutron fluid region. Here the density could reach somewhere between $3\times 10^{17}$ and $\sim 10^{18}$ kg/m$^3$ before something may (or may not) cause a phase change (hyperons, quark matter etc.). The pressure at the upper end of this range would approach that which I gave above. A more formal calculation could calculate the equilibrium composition of an n,p,e gas and then use the appropriate formula for an ideal Fermi gas to estimate the total pressure. Such a calculation, at a range of densities, also gives you an estimate of the adiabatic index. I find that the pressure of an ideal n,p,e gas is $7\times 10^{32}$ Pa at $3\times 10^{17}$ kg/m$^3$ (completely dominanted by non-relativistic neutrons), rising to $5\times 10^{33}$ Pa at $10^{18}$ kg/m$^3$. The adiabatic index $\gamma$, (where I define $\gamma$ in terms of $P \propto \rho^{\gamma}$), is approximately 1.6 over this range.
You have to realise though that this is an unrealistic equation of state at these densities. As the neutrons get closer together ($\leq 10^{-15}$ m) they must interact and this hardens the equation of state such that $\gamma \geq 2$, otherwise neutron stars of mass up to $2M_{\odot}$ couldn't exist! Indeed, one definition of the "inner core" might be where this hardening happens.
A recent paper by Hebeler et al. (2013) includes the following plot (in cgs units) which reviews the many possible equations of state in the density range I have discussed above. For reasons discussed in that paper, the authors favour the AP3 or AP4 equations of state which give higher pressures tan the simple n,p,e approximation. The slope of these graphs gives the the value of $\gamma \simeq 3.3$ for AP3 and AP4 for $\rho > 5 \times 10^{17}$ kg/m$^3$.
• @umarkhan It will give you a result, but the outer parts of a neutron star would not be well characterised by such an equation of state. $\gamma$ is smaller in the n,p,e fluid below nuclear densities and in the crust. Nov 8, 2016 at 9:41


text:
# (solved)Question 2.34 of NCERT Class XI Chemistry Chapter 2

Calculate the energy required for the process He+ (g) → He2+ (g) + e- The ionization energy for the H atom in the ground state is 2.18 × 10-18 J atom-1
(Rev. 19-Dec-2022)

## Categories | About Hoven's Blog

,

### Question 2.34NCERT Class XI Chemistry

Calculate the energy required for the process He+ (g) → He2+ (g) + e- The ionization energy for the H atom in the ground state is 2.18 × 10-18 J atom-1

### Solution in Detail(video solution below this)

by Bohr's theory, IE is the energy
required to remove an electron from
the n-th shell of a H-like atom of
atomic number Z:

$\displaystyle E_n = +R_H\bigg(\frac{Z^2}{n^2}\bigg)$ J/atom

∴ IE of hydrogen for n = 1 and Z = 1 is

$\displaystyle E_1 = R_H \bigg(\frac{1^2}{1^2}\bigg) = R_H \text{ (J/atom)}$

but given $\displaystyle E_1 = 2.18 \times 10^{-18}\text{ J/atom}$

$\displaystyle \therefore R_H = 2.18 \times 10^{-18}\text{ . . . (1)}$

He+ is hydrogen-like with 1 electron

hence, IE for He+ (Z = 2, n = 1)

$\displaystyle \therefore \text{ IE } = +R_H\bigg(\frac{2^2}{1^2}\bigg) = R_H \times 4$

$\displaystyle = 2.18\times10^{-18} \times 4 \:\bigg[\text{R}_H \text{ by (1)}\bigg]$

$\displaystyle = 8.72 \times 10^{-18} \text{ J/atom }\underline{Ans}$

### Video Explanation

Please watch this youtube video for a quick explanation of the solution:


text:
Intermediate Algebra (12th Edition)

$(a-4)(2a+1)$
$\bf{\text{Solution Outline:}}$ To factor the given expression, $2a^2-7a-4 ,$ use the factoring of trinomials in the form $ax^2+bx+c.$ $\bf{\text{Solution Details:}}$ In the trinomial expression above, $a= 2 ,b= -7 ,\text{ and } c= -4 .$ Using the factoring of trinomials in the form $ax^2+bx+c,$ the two numbers whose product is $ac= 2(-4)=-8$ and whose sum is $b$ are $\left\{ -8,1 \right\}.$ Using these two numbers to decompose the middle term results to \begin{array}{l}\require{cancel} 2a^2-8a+1a-4 .\end{array} Using factoring by grouping, the expression above is equivalent to \begin{array}{l}\require{cancel} (2a^2-8a)+(1a-4) \\\\= 2a(a-4)+(a-4) \\\\= (a-4)(2a+1) .\end{array}


text:
# “Understanding Machine Learning: From Theory to Algorithms” the universal approximation theorem

I'm readying on "Understanding Machine Learning: From Theory to Algorithms" the Universal approximation theorem:

..."Networks are universal approximators. That is, for every fixed precision parameter, $$\epsilon >0$$, and every Lipschitz function $$f : [−1; 1]^{n} \rightarrow [−1; 1]$$, it is possible to construct a network such that for every input $$\textbf{x} \in [−1; 1]^{n}$$, the network outputs a number between $$f(x) − \epsilon$$ and $$f(x) + \epsilon$$".

It seems to me that the function $$f : [−1; 1]^{n} \rightarrow [−1; 1]$$ is Boolean.

I think that the n-dimensional unit hypercube $${\displaystyle [-1,1]^{n}}$$ can be replaced with a compact set of $$\mathbb{R} ^ {n}$$ but the codomain makes me puzzled.

I was expecting a function: $$f : \mathbb { R } ^ { n } \rightarrow \mathbb { R }$$

• I think the question is best suited for CrossValidated. – Shubham Panchal Jul 11 '19 at 3:16
• Did you mean to use semicolons in the intervals? The book you're referring to uses commas throughout. Are you using semicolons to indicate the discrete set (hence "is Boolean"? I'm not familiar with that notation... – Ben Reiniger Jul 11 '19 at 15:14

Yes, it is applicable to a Boolean function and you can modify the result to compact set.

A result where the codomian is $$[-1, 1]$$ would cover regression and not restricting ourself to boolean.

You might or might not have a result for $$f: \mathbb{R}^n \to \mathbb{R}$$. A theorem states what conditions it requires in order to hold, for conditions that are not promised by a theorem, you have to check carefully whether the same conclusion can still apply. In particular, if the proof uses compactness, then you can't replace $$[-1,1]$$ with $$\mathbb{R}$$.

There are different versions of the theorem; see wikipedia's for example.

Some counterexamples will help focus the boundaries of what is and what is not possible. I won't be too formal with this, but I think it could be made so with some effort. Let's assume we're using the ReLU activation function, so that everything in the end is piecewise-linear and continuous. With a fixed network architecture, the number of hinge points in the final function is bounded (in a network with just one hidden layer, it is at most twice the number of hidden neurons). So, $$f: \mathbb{R}\to[-1,1],\ f(x)=\sin(x)$$ cannot be approximated: beyond all hinge points we have a single linear function, which cannot wiggle like $$\sin$$. You'll have the same problem with an open domain, using something akin to $$\sin(1/x)$$. So we need a compact domain.

You also obviously need continuity if you want approximation in the sense that every input is close to the right output: a jump discontinuity can be approximated with a very high slope, but only in the sense that the measure of the poorly-approximated points can be made as small (but positive) as desired.

And now, the continuous image of a compact set is compact, so the range is necessarily compact already. (So here you could safely expand the theorem to have codomain $$\mathbb{R}$$.)


text:
# Number guessing game with expensive answers

I play the number guessing game with positive integers and a known constant upper bound. Every time I make a guess about the number I have to pay 1 dollar but if my partner answers 'Yes' I have to pay 9 more.

What is the best strategy to minimize my average cost of the game (other than not playing)?

Does the strategy change if I have to pay more/less for the good guesses?

• Exchange rôles with your partner! :-) Seriously, does the game continue until you guess the number? – Brian M. Scott Sep 9 '13 at 22:56
• If you guess wrong, are you given any information, like higher or lower? Otherwise, do you have any information about the probability distribution? If not, you will have to guess half the numbers on average. – Ross Millikan Sep 9 '13 at 23:20
• I am given an answer to every guess (and my partner doesn't lie), the trick is that my cost depends on the answer. Example: Is the number greater than 10? Yes. => I pay 10. Is the number greater than 12? No => I pay 1. The game goes on until I guess the exact number. If I have info about the distribution I should simply guess for the more probable values first right? – buherator Sep 10 '13 at 5:59

Ordering doesn't matter here, only the size of the search space. Since every question is yes/no, it splits the space into two parts, and the only parameter that matters for each question is how to size those parts. This allows us to compute the expected cost recursively. Define $T(n)$ to be the average cost of asking enough questions to determine the number from a search space of $n$. By definition $T(1)=0$ (technically the game doesn't end until you have actually guessed the number, but this adds only $O(1)$ to the cost so I'll play fast and loose with this distinction).

For $n \ge 2$, the expected cost of asking "is the number in a given subset of size $k$?" and continuing to play the game on the remaining search space is $1 + \frac{k}{n} (9 + T(k)) + \frac{n-k}{n} T(n-k)$. So we optimize this quantity over all $0 < k < n$:

$$T(n) = \min_{0<k<n} 1 + \frac{k}{n} (9 + T(k)) + \frac{n-k}{n} T(n-k) = 1 + \frac1n \min_{0<k<n} (9k + kT(k) + (n-k)T(n-k)).$$

This gives an optimal value of $T(100) = 26.5$, substantially better than CJ from buherator's answer (assuming we agree on the $O(1)$ convention). I computed the optimal costs and moves up to $n=10000$.

An optimal strategy for $n=10000$ costs about $52.1$, and the first question is to ask "is the number in the range $[1,1607]$?" This has a $0.1607$ chance of costing $10$ plus the cost of solving the size $1607$ problem (which is $\approx 41.956$), and a $0.8393$ chance of costing $1$ plus the cost of solving the size $8393$ problem (which is $\approx 51.127$).

This example suggests that the optimal cutoff is approximately where the two cases balance out in cost, just as in the symmetric case. If we make this assumption then we can do a pretty good job of approximating the optimal strategy. The value of $T(n)$ is a priori somewhere between $\log_2 n$ and $10 \log_2 n$: we make the ansatz that it is asymptotic to $c \log n$ for some constant $c$. Then the optimal cutoff occurs where

$$9 + c \log k = c \log(n-k) \implies \frac{k}{n} = \frac{1}{1 + e^{9/c}}.$$

On the other hand (still under the key assumption) the cost is also equal to the number of iterations that descend into the large half, so $c = -1/\log(1-\epsilon)$. This gives two equations in $\epsilon$ and $c$:

$$\epsilon = \frac{1}{1 + e^{9/c}}, c = -1/\log(1-\epsilon),$$ which has the numerical solution $\epsilon \approx 0.164922, c \approx 5.548537$. So I would conjecture that the optimal strategy is close to "is the number in the first 16.5% of the current range?" The actual ratio appears to be rather non-monotonic, varying up and down between $0.14$ and $0.17$ as $n$ increases, even for large $n$. This could be due to rounding errors in my computation (we could be comparing several $k$ values whose scores are really close), but if not then it suggests that the deviation between the optimal $k$ and $\epsilon n$ is more than just rounding to an integer value of $k$.

• Thank you for your detailed answer! I'm trying to implement some experiments to test it (I'm not smart enough for a formal proof..), do I understand correctly that if the penalty cost changes only the $\epsilon$ part changes (in e's exponent) in the final equation system? – buherator Oct 6 '16 at 12:22
• @buherator Yeah I think the relative penalty appears as the $9$ throughout. Of course the final pair of equations is only conjectural. – Erick Wong Oct 6 '16 at 12:34
• Your solution seems to give convincing results with different parameters too: github.com/v-p-b/campzer0-numberguess . I'm accepting your answer, thank you! – buherator Oct 11 '16 at 15:38
• @buherator By the way, the first formula should be exactly optimal. The conjecture is that for large $n$ it converges to that value of $\epsilon$. It would be interesting to see how the two strategies compare to each other. – Erick Wong Oct 11 '16 at 18:10
• (The rounding issue I discuss in the recent edit only affects the computation of the exact optimal value of $k$... it should have only a minimal effect on the simulated or theoretical cost). – Erick Wong Oct 11 '16 at 18:42

So after all, we put together a little sample script for simulation:

https://gitorious.org/campzer0/numberguessing/

We implemented the naive linear approach and two modified version of binary search:

• synalgo performs simple binary search and tries to estimate the expected remaining cost. If the expected cost of moving on with binary search is higher than the expected cost of the linear search it falls back to linear
• cj performs a modified binary search where the search space is split in the ration of the good answer/bad answer costs - in the case described in the questions the algorithm always asks if the number is in the lower 1/10 of the current interval. This algorithm also falls back to linear search when the 1/10th of the interval becomes less than 1

After 10000 runs with an upper bound of 100 the results look like this:

naive average of 10000 runs: 59.8012
synalgo average of 10000 runs: 32.9129
cj average of 10000 runs: 30.9006

CJ is usually ~10% better than synalgo and they both present good results compared to the naive approach.

These algorithms are of course not proven (close-to) optimal, but provide acceptable efficiency for my actual problem. Any further optimization proposals are welcome though!


text:
# reMarkable 2 review

I recently got a reMarkable 2. I’ve had it for one week. This is my review of it so far.

## TLDR

• It’s very very good and I’m happy I bought one. If your relationship with working “on things” is like mine, I recommend it.
• It’s expensive and you can get an iPad + accessories + Apple Pencil for the same price, which may be better
• There are some annoyances

Would I buy the reMarkable again? That is, if it was struck by a meteor tomorrow, would I buy a new one? Yes (but see my notes on other options). In my opinion, these devices are now good enough that it’s worth paying a premium to get the “paper experience”. I would not recommend saving up for the reMarkable - that is, if your disposable income is small and you’d need to make an effort over a time period to find the money, it’s probably not worth it. If you own a laptop and are considering buying a tablet, which you anticipate using largely for reading and writing, the reMarkable is recommendable. If you’re in the market for a tablet to use also for laptop tasks, it’s not up to the job.

## Paper and Me

My fiancé is one of those stationery nerds. She own several different fountain pens and a metric boatload of ballpoint pens, gel pens, highlighters, markers, and what have you. She buys a specific brand of Japanese calendar (Hobonichi) each year. She has had the same color-coding system for her notes since we started university 1. She meticulously pagemarks and sorts all her notes into plastic sleeves and binders.

I am not like my fiancé. My workflow for most of my first five years at university consisted of the following:

• Keep a stack of papers and ballpoint pens in my backpack (along with the books I’m using and my laptop and so on).
• When you sit down to work somewhere, get laptop, stack of papers, and a ballpoint pen out.
• When you need to think, locate a blank spot on a piece of paper in front of you and write there.
• If there’s no space, go through the papers until you find some space.
• When you stop working, gather all the papers into your backpack.
• Periodically go through your papers and throw out any that have writing all over them
• Periodically add a fresh dump of blank papers to your backpack.

The two points I want to make here is:

• I am a very disorganized person
• I fucking love paper

Paper is awesome. I am long paper. Whoever you are, you should be using more paper. There’s absolutely no reason to strain your head to the breaking point, keeping a bunch of things inside it, instead of putting them on paper. Actually, if you take one recommendation from this review, just make sure to keep paper and pen wherever you work and use it all the time. Go to these links to buy some right now, they’re cheap (but you probably have this stuff already - use it!). See also Paper Trauma. I also invite you to consider buying a big notebook, but I haven’t tried that yet myself.

But paper has some limitations:

• You sometimes run out
• If you want to write your notes about something you’re reading on paper, the best place is the thing you’re reading - but this requires that it be on paper (instead of a computer). This means you have to buy a book (these have a tendency to pile up rapidly, if it even is a book you can buy) or print it out (this tends to be messy, and you might not have a printer. Back at the university of Copenhagen, I used to print out a lot of stuff and spiral-back it, but I can’t do this any more).
• Sometimes you lose your old paper where you wrote something important (usually, because you didn’t know at the time it would be important).
• It’s easy to make a mess by creating a huge stack of paper. If you’re like me, “clean up your big pile of papers” can sometimes become a bit of an “ugh” task that you can’t get around to and which sucks up a lot of energy.

Enter the reMarkable 2. A “next-generation paper tablet”, it promises you that you can “Replace your notebooks and printed documents with the only tablet that feels like paper”. Naturally, I was intrigued.

## General experience and thoughts

The writing experience is very very good. It doesn’t quite feel like paper - it’s not that close - but apart from actually writing on paper it’s the best I’ve tried (having, admittedly, not tried that many). It’s leagues ahead of writing on a “normal” tablet using a stylus.

Drawing on the reMarkable has the following advantages over paper:

• You can erase things, perfectly
• You can undo/redo strokes
• You can zoom
• You can select and move things
• There’s layers, like in image editing software.

The surface feels a bit too smooth - it doesn’t have the same friction as dragging a pen over paper. It’s still better than a normal screen (like a iPad) though. You can tell they added some texture to make this more lifelike, but it’s not the same.

Not surprisingly, it’s a bit awkward to hold the tablet with one hand and write with the other. The holding hand doesn’t provide quite enough support. You really need to support it with a table or your knees or something.

Probably the biggest thing I would like is a bigger screen.

The “OS”, as it is, is basically just a folder structure where you can arrange pdfs, EPUBs (I haven’t tried this), and “notebooks”, which are basically collections of pages of your drawing/writing/scribbling. You can also add things to “favorites” and sort by “type” (notebook/pdf/EPUB). It’s not amazing, especially since it’s a bit sluggish to use, but it does the job.

You can get stuff out of and into the reMarkable in a few ways:

• If you create an account in their server and connect it, it automatically syncs to their cloud, and you can get at this data from a desktop app and a mobile app. These are also how you get stuff into it. The syncing here is pretty slow, so I don’t use this a lot.

• You can send a mail from the tablet with a collection of pages - either annotated pdf/epub or raw notebook. This is what I do, and it works great!

• There’s also a browser extension which turns a website into an EPUB and sends it to the tablet. It works okay, about the same as something like Pocket or Instapaper.

• You can also have the reMarkable convert you handwriting into text and send it as an email. This is actually surprisingly (to me) good, but probably not quite good enough to send a mail without proofreading (you can check and edit the message on the tablet before sending). The first line in the following picture was converted without error, while the second became “The Quick brown tote jumps over the lazy chef”.

I haven’t tested this extensively, and I expect it to do worse on other languages than English (you can switch the language used in the setting menu).

### The physical feel of things

The pen weighs 19g, and feels about like an ordinary ballpoint pen/fineliner in the hand. I would have liked a bit more weight, but it’s not a big deal.

The tablet is pretty light. Holding it in one hand is a bit heavy for me - not TOO heavy to use like that, but I probably won’t read like that for, say, multiple hours. This is not a big deal though.

The “Book Folio” cover is fine. It’s fairly basic, certainly not “worth” 100 EUR - if there was a market for accesories, you could certainly buy a better one much cheaper. It essentially consists of two pieces of sturdy cardboard (I think), covered on the outside by light grey fabric and on the inside by something that feels like felt. The tablet attached to the backside by magnets - these are not that strong, but they do the job. The pen simply attached to the tablet itself by magnets - again, these do the job, but I would have liked them to be a bit stronger. I am slightly worried about the pen falling off in a bag full of mess and getting into trouble. The front piece also attaches by (weaker) magnets when closed, so it doesn’t open by itself.

The basic things I’ve done with the reMarkable so far breaks down into roughly three categories:

• During video meetings, I’ve scratched notes on the tablet
• I’ve put things on the tablet to read, and scribbled in the margins while I did so
• I’ve used the blank pages (notebooks) to work on.

### Meeting notes

The device doesn’t really shine here. The only advantage it has over writing on the computer is that I can keep the window with people’s faces visible while writing, which I could do anyway if I had a second monitor. There’s something nice about paper writing - maybe the flexibility of not being limited to one line after the other? Maybe something physical? That this captures but writing on the computer doesn’t, though. But this is a small advantage.

This is the real killer app for me. As noted above, you can put your pdfs on the reMarkable and read them. And while reading them, you can scribble on them. I am one of those people who likes to constantly doodle on things, so the eraser really comes in handy here - I can make a scribble, then erase it, so the margins don’t get covered in random noise. This is an advantage I hadn’t thought of.

After scribbling and highlighting and so on, I usually send annotated pages to my email. Then they can become part of my diginal notes.

I think the fact that this doesn’t browse the internet is probably a huge advantage, because it adds a bit of friction to checking twitter, making me distract myself less. Just having a pdf reader that can’t check twitter might contribute a huge amount of value.

### Working on blank paper

As mentioned above, using blank paper to think on is really important to me. On the reMarkable, there’s a distinguished “quick pages” notebook with a permanent shortcut to it on the top. This gives you an infinite supply of blank pages to work on. I also use a few other notebooks to collect thoughts about various projects, but this is not very well-organized. The biggest drawback here is that the pages are not very large. I think the screen is slightly broader than A6 format (checked by holding up a folded sheet of A4). This is not optimal - I would like to be able to draw bigger things, to keep more thoughts on the page at once. The jury is still out, but I can’t imagine the reMarkable will permanently replace paper in this domain for me - at best, it will supplement it.

## Notes on price & comparison to other products

The reMarkable 2 device itself costs 400 EUR (473 USD at current rates, but prices may differ by region). This price point is severely misleading, since it’s more or less useless unless you buy a pen (you can read on it, of course, but if that’s all you want to do I would definitely consider it overpriced).

The basic “Marker” costs 60 EUR, and the “Marker plus” costs 100 EUR. The only difference is that the plus weighs slightly more and that you can use the rear end to erase. It’s probably best to just consider the base price of the device 460 EUR and ask whether you want to pay 40 EUR for the eraser.

Lastly, you can purchase a cover, or “folio”. Here your options are the basic folio for 80 EUR, which is simply a sleeve that you can slide the reMarkable into, and the “Book Folio”, which is essentially a folder, which magnetically attaches to the back of the reMarkable and wraps around it. This costs 100 EUR for the basic “polymer weave”, and 150 EUR for the black or brown “Premium Leather” version. These are “overpriced” in the sense that, if there was a market for reMarkable accessories, you would certainly be able to buy much better protective covers for much less. But there isn’t, so you’re stuck with these, and the question is whether they’re worth the money.

You can also buy a reMarkable 1 with basic Folio and (non-eraser) Marker for 349 EUR.

When I bought my reMarkable, I paid a reduced price because I was buying it early. I bought the Marker Pro and the Polymer Weave Book Folio, and ended up paying 460 EUR.

So how much does the reMarkable 2 cost? Based on my experience so far, the built-in eraser is super convenient and well worth the money. You probably want some sort of protective cover, to prevent the screen from getting dings and scratches, but I have to stress that the Folios really seem to be nothing special. If I was buying it again I would probably buy the Book Folio again, just because 20 EUR isn’t a lot and it seems a bit incovenient to be taking it in and out of the sleeve all the time. On that basis, it costs 600 EUR, or ~710 USD.

That’s a lot. For that money you could buy an iPad (329 USD), Apple Pencil ($129), and a smart cover ($50), and still have money left over. You could use that money to buy extra memory for the iPad, cellular web access (not a killer app for me, but the reMarkable lacks it). You could even buy an iPad Air, if you maybe skimped out on the cover. I didn’t look into the Android options, but presumably those are even cheaper. So if you’re looking for a “general work tablet”, you should probably skip the reMarkable. If you’re looking for the paper feeling, you could consider Paperlike, for just 40 USD (which I have not tried and cannot opine on).

reMarkable really has nothing in the way of tablet features. You can’t install any apps (actually, there’s a small community of hackers, but nobody is seriously developing third-party software). You wouldn’t want to install any of the apps you have on other devices even if you could, because the refresh speed is too low for most “normal” software to really be usable. You can’t browse the internet.

What you can do is basically:

• Write and draw with a stylus like you were drawing on blank paper
• Read pdfs and ebooks, and write on them.

If stuff like that is such an important part of your workflow that getting a significantly better experience than you would doing the same thing on an iPad is worth giving up all of the iPad’s other features, you should consider the reMarkable. Personally, this is basically the only reason I would really want to own a tablet anyway. In fact, the fact that I can’t browse the internet on it is probably a strong feature in its favor, since it makes it harder to distract yourself with twitter.

Another comparison is with the line of products from Boox. The most direct comparison is with the Note Air, which has a similar form factor and price point. There is also Supernote. I have no experience with any of these. It seems that the Boox Note Air is generally more feature-rich than the reMarkable - it actually runs Android, and is thus closer to a “proper” tablet. This review may provide some information, and recommends the Note Air.

I have to say these products are really tempting me, especially the ones with larger screens. If I was buying a new “paper tablet” today, I would seriously consider them.

1. for the record, her opinion of my reMarkable was that it seemed awesome, but it was a nonstarter for her, because it didn’t do colors ↩


text:
A variational approach to stable principal component pursuit

# A variational approach to stable principal component pursuit

## Abstract

We introduce a new convex formulation for stable principal component pursuit (SPCP) to decompose noisy signals into low-rank and sparse representations. For numerical solutions of our SPCP formulation, we first develop a convex variational framework and then accelerate it with quasi-Newton methods. We show, via synthetic and real data experiments, that our approach offers advantages over the classical SPCP formulations in scalability and practical parameter selection.

## 1 Introduction

Linear superposition is a useful model for many applications, including nonlinear mixing problems. Surprisingly, we can perfectly distinguish multiple elements in a given signal using convex optimization as long as they are concise and look sufficiently different from one another. Popular examples include robust principal component analysis (RPCA) where we decompose a signal into low rank and sparse components and stable principal component pursuit (SPCP), where we also seek an explicit noise component within the RPCA decomposition. Applications include alignment of occluded images (Peng et al., 2012), scene triangulation (Zhang et al., 2011), model selection (Chandrasekaran et al., 2012), face recognition, and document indexing (Candès et al., 2011).

The SPCP formulation can be mathematically stated as follows. Given a noisy matrix , we decompose it as a sum of a low-rank matrix and a sparse matrix via the following convex program

 Missing dimension or its units for \kern ($\text{SPCP}_{\text{sum}}$) subject to∥L+S−Y∥F≤ε,

where the 1-norm and nuclear norm are given by where is the vector of singular values of . In (($\text{SPCP}_{\text{sum}}$)), the parameter controls the relative importance of the low-rank term vs. the sparse term , and the parameter accounts for the unknown perturbations in the data not explained by and .

When , (($\text{SPCP}_{\text{sum}}$)) is the “robust PCA” problem as analyzed by Chandrasekaran et al. (2009); Candès et al. (2011), and it has perfect recovery guarantees under stylized incoherence assumptions. There is even theoretical guidance for selecting a minimax optimal regularization parameter (Candès et al., 2011). Unfortunately, many practical problems only approximately satisfy the idealized assumptions, and hence, we typically tune RPCA via cross-validation techniques. SPCP further complicates the practical tuning due to the additional parameter .

To cope with practical tuning issues of SPCP, we propose the following new variant called “max-SPCP”:

 minimizeL,Smax(|||L|||∗,λmax∥S∥1) ($\text{SPCP}_{\text{max}}$) subject to∥L+S−Y∥F≤ε,

where acts similar to . Our work shows that this new formulation offers both modeling and computational advantages over (($\text{SPCP}_{\text{sum}}$)).

Cross-validation with (($\text{SPCP}_{\text{max}}$)) to estimate is significantly easier than estimating in (($\text{SPCP}_{\text{sum}}$)). For example, given an oracle that provides an ideal separation , we can use in both cases. However, while we can estimate , it is not clear how to choose from data. Such cross validation can be performed on a similar dataset, or it could be obtained from a probabilistic model.

Our convex approach for solving (($\text{SPCP}_{\text{sum}}$)) generalizes to other source separation problems (Baldassarre et al., 2013) beyond SPCP. Both (($\text{SPCP}_{\text{max}}$)) and (($\text{SPCP}_{\text{sum}}$)) are challenging to solve when the dimensions are large. We show in this paper that these problems can be solved more efficiently by solving a few (typically 5 to 10) subproblems of a different functional form. While the efficiency of the solution algorithms for (($\text{SPCP}_{\text{sum}}$)) relies heavily on the efficiency of the 1-norm and nuclear norm projections, the efficiency of our solution algorithm (($\text{SPCP}_{\text{max}}$)) is preserved for arbitrary norms. Moreover, (($\text{SPCP}_{\text{max}}$)) allows a faster algorithm in the standard case, discussed in Section 6.

## 2 A Primer on Spcp

The theoretical and algorithmic research on SPCP formulations (and source separation in general) is rapidly evolving. Hence, it is important to set the stage first in terms of the available formulations to highlight our contributions.

To this end, we illustrate (($\text{SPCP}_{\text{sum}}$)) and (($\text{SPCP}_{\text{max}}$)) via different convex formulations. Flipping the objective and the constraints in (($\text{SPCP}_{\text{max}}$)) and (($\text{SPCP}_{\text{sum}}$)), we obtain the following convex programs

 minimizeL,S12∥L+S−Y∥2F ($\text{flip-SPCP}_{\text{sum}}$) s.t.|||L|||∗+λsum∥S∥1≤τsum
 minimizeL,S12∥L+S−Y∥2F ($\text{flip-SPCP}_{\text{max}}$) Missing or unrecognized delimiter for \left
###### Remark 2.1.

The solutions of (($\text{flip-SPCP}_{\text{sum}}$)) and (($\text{flip-SPCP}_{\text{max}}$)) are related to the solutions of (($\text{SPCP}_{\text{sum}}$)) and (($\text{SPCP}_{\text{max}}$)) via the Pareto frontier by Aravkin et al. (2013a, Theorem 2.1). If the constraint is tight at the solution, then there exist corresponding parameters and , for which the optimal value of (($\text{flip-SPCP}_{\text{sum}}$)) and (($\text{flip-SPCP}_{\text{max}}$)) is , and the corresponding optimal solutions and are also optimal for (($\text{SPCP}_{\text{sum}}$)) and (($\text{SPCP}_{\text{max}}$)).

For completeness, we also include the Lagrangian formulation, which is covered by our new algorithm:

 minimizeL,SλL|||L|||∗+λS∥S∥1+12∥L+S−Y∥2F (Lag-SPCP)

Problems (($\text{flip-SPCP}_{\text{max}}$)) and (($\text{flip-SPCP}_{\text{sum}}$)) can be solved using projected gradient and accelerated gradient methods. The disadvantage of some of these formulations is that it may not be clear how to tune the parameters. Surprisingly, an algorithm we propose in this paper can solve (($\text{SPCP}_{\text{max}}$)) and (($\text{SPCP}_{\text{sum}}$)) using a sequence of flipped problems that specifically exploits the structured relationship cited in Remark 2.1. In practice, we will see that better tuning also leads to faster algorithms, e.g., fixing ahead of time to an estimated ‘noise floor’ greatly reduces the amount of required computation if parameters are to be selected via cross-validation.

Finally, we note that in some cases, it is useful to change the term to where is a linear operator. For example, let be a subset of the indices of a matrix. We may only observe restricted to these entries, denoted , in which case we choose . Most existing RPCA/SPCP algorithms adapt to the case but this is due to the strong properties of the projection operator . The advantage of our approach is that it seamlessly handles arbitrary linear operators . In fact, it also generalizes to smooth misfit penalties, that are more robust than the Frobenius norm, including the Huber loss. Our results also generalize to some other penalties on besides the 1-norm.

The paper proceeds as follows. In Section 3, we describe previous work and algorithms for SPCP and RPCA. In Section 4, we cast the relationships between pairs of problems (($\text{flip-SPCP}_{\text{sum}}$)), (($\text{SPCP}_{\text{sum}}$)) and (($\text{flip-SPCP}_{\text{max}}$)), (($\text{SPCP}_{\text{max}}$)) into a general variational framework, and highlight the product-space regularization structure that enables us solve the formulations of interest using corresponding flipped problems. We discuss computationally efficient projections as optimization workhorses in Section 5, and develop new accelerated projected quasi-Newton methods for the flipped and Lagrangian formulations in Section 6. Finally, we demonstrate the efficacy of the new solvers and the overall formulation on synthetic problems and a real cloud removal example in Section 7, and follow with conclusions in Section 8.

## 3 Prior Art

While problem (($\text{SPCP}_{\text{sum}}$)) with has several solvers (e.g., it can be solved by applying the widely known Alternating Directions Method of Multipliers (ADMM)/Douglas-Rachford method (Combettes & Pesquet, 2007)), the formulation assumes the data are noise free. Unfortunately, the presence of noise we consider in this paper introduces a third term in the ADMM framework, where the algorithm is shown to be non-convergent (Chen et al., 2013). Interestingly, there are only a handful of methods that can handle this case. Those using smoothing techniques no longer promote exactly sparse and/or exactly low-rank solutions (Aybat et al., 2013). Those using dual decomposition techniques require high iteration counts. Because each step requires a partial singular value decomposition (SVD) of a large matrix, it is critical that the methods only take a few iterations.

As a rough comparison, we start with related solvers that solve (($\text{SPCP}_{\text{sum}}$)) for . Wright et al. (2009a) solves an instance of (($\text{SPCP}_{\text{sum}}$)) with and a system in hours. By switching to the ((Lag-SPCP)) formulation, Ganesh et al. (2009) uses the accelerated proximal gradient method (Beck & Teboulle, 2009) to solve a matrix in under one hour. This is improved further in Lin et al. (2010) which again solves (($\text{SPCP}_{\text{sum}}$)) with using the augmented Lagrangian and ADMM methods and solves a system in about a minute. As a prelude to our results, our method can solve some systems of this size in about seconds (c.f., Fig. 1).

In the case of (($\text{SPCP}_{\text{sum}}$)) with , Tao & Yuan (2011) propose the alternating splitting augmented Lagrangian method (ASALM), which exploits separability of the objective in the splitting scheme, and can solve a system in about five minutes.

The partially smooth proximal gradient (PSPG) approach of Aybat et al. (2013) smooths just the nuclear norm term and then applies the well-known FISTA algorithm (Beck & Teboulle, 2009). Aybat et al. (2013) show that the proximity step can be solved efficiently in closed-form, and the dominant cost at every iteration is that of the partial SVD. They include some examples on video, lopsided matrices: or so, in about 1 minute). solving formulations in under half a minute.

The nonsmooth adaptive Lagrangian (NSA) algorithm of Aybat & Iyengar (2014) is a variant of the ADMM for (($\text{SPCP}_{\text{sum}}$)), and makes use of the insight of Aybat et al. (2013). The ADMM variant is interesting in that it splits the variable , rather than the sum or residual . Their experiments solve a 1500 1500 synthetic problems in between 16 and 50 seconds (depending on accuracy) .

Shen et al. (2014) develop a method exploiting low-rank matrix factorization scheme, maintaining . This technique has also been effectively used in practice for matrix completion (Lee et al., 2010; Recht & Ré, 2011; Aravkin et al., 2013b), but lacks a full convergence theory in either context. The method of (Shen et al., 2014) was an order of magnitude faster than ASALM, but encountered difficulties in some experiments where the sparse component dominated the low rank component in some sense. Mansour & Vetro (2014) attack the ($\text{SPCP}_{\text{sum}}$) formulation using a factorized approach, together with alternating solves between and . Non-convex techniques also include hard thresholding approaches, e.g. the approach of Kyrillidis & Cevher (2014). While the factorization technique may potentially speed up some of the methods presented here, we leave this to future work, and only work with convex formulations.

## 4 Variational Framework

Both of the formulations of interest (($\text{SPCP}_{\text{sum}}$)) and (($\text{SPCP}_{\text{max}}$)) can be written as follows:

 minϕ(L,S)s.t.ρ(L+S−Y)≤ε. (4.1)

Classic formulations assume to be the Frobenius norm; however, this restriction is not necessary, and we consider to be smooth and convex. In particular, can be taken to be the robust Huber penalty (Huber, 2004). Even more importantly, this formulation allows pre-composition of a smooth convex penalty with an arbitrary linear operator , which extends the proposed approach to a much more general class of problems. Note that a simple operator is already embedded in both formulations of interest:

 L+S=[II][LS]. (4.2)

Projection onto a set of observed indices is also a simple linear operator that can be included in . Operators may include different transforms (e.g., Fourier) applied to either or .

The main formulations of interest differ only in the functional . For (($\text{SPCP}_{\text{sum}}$)), we have

 ϕsum(L,S)=|||L|||∗+λsum∥S∥1,

while for (($\text{SPCP}_{\text{max}}$)),

 ϕmax(L,S)=max(|||L|||∗,λmax∥S∥1).

The problem class (4.1) falls into the class of problems studied by van den Berg & Friedlander (2008, 2011) for and by Aravkin et al. (2013a) for arbitrary convex . Making use of this framework, we can define a value function

 v(τ)=minL,Sρ(A(L,S)−Y)s.t. ϕ(L,S)≤τ, (4.3)

and use Newton’s method to find a solution to . The approach is agnostic to the linear operator (it can be of the simple form (4.2); include restriction in the missing data case, etc.).

For both formulations of interest, is a norm defined on a product space , since we can write

 ϕsum(L,S) =∥∥∥|||L|||∗λsum∥S∥1∥∥∥1, (4.4) ϕmax(L,S) =∥∥∥|||L|||∗λmax∥S∥1∥∥∥∞. (4.5)

In particular, both and are gauges. For a convex set containing the origin, the gauge is defined by

 γ(x∣C)=infλ{λ:x∈λC}. (4.6)

For any norm , the set defining it as a gauge is simply the unit ball . We introduce gauges for two reasons. First, they are more general (a gauge is a norm only if is bounded with nonempty interior and symmetric about the origin). For example, gauges trivially allow inclusion of non-negativity constraints. Second, definition (4.6) and the explicit set simplify the exposition of the following results.

In order to implement Newton’s method for (4.3), the optimization problem to evaluate must be solved (fully or approximately) to obtain . Then the parameter for the next (4.3) problem is updated via

 τk+1=τk−v(τ)−τv′(τ). (4.7)

Given , can be written in closed form using (Aravkin et al., 2013a, Theorem 5.2), which simplifies to

 v′(τ)=−ϕ∘(AT∇ρ(A(¯L,¯¯¯¯S)−Y)), (4.8)

with denoting the polar gauge to . The polar gauge is precisely , with

 C∘={v:⟨v,x⟩≤1∀x∈C}. (4.9)

In the simplest case, where is given by (4.2), and is the least squares penalty, the formula (4.8) becomes

 v′(τ)=−ϕ∘([¯¯¯¯L+¯¯¯¯S−Y¯¯¯¯L+¯¯¯¯S−Y]).

The main computational challenge in the approach outlined in (4.3)-(4.8) is to design a fast solver to evaluate . Section 6 does just this.

The key to RPCA is that the regularization functional is a gauge over the product space used to decompose into summands and . This makes it straightforward to compute polar results for both and .

###### Theorem 4.1 (Max-Sum Duality for Gauges on Product Spaces).

Let and be gauges on and , and consider the function

 g(x,y)=max{γ1(x),γ2(y)}.

Then is a gauge, and its polar is given by

 g∘(z1,z2)=γ∘1(z1)+γ∘2(z2).
###### Proof.

Let and denote the canonical sets corresponding to gauges and . It immediately follows that is a gauge for the set , since

 inf{λ≥0|(x,y)∈λC} =inf{λ|x∈λC1 and y∈λC2} =max{γ1(x),γ2(y)}.

By (Rockafellar, 1970, Corollary 15.1.2), the polar of the gauge of is the support function of , which is given by

 supx∈C1,y∈C2⟨(x,y),(z1,z2)⟩ =supx∈C1⟨x,z1⟩+supy∈C2⟨y,z2⟩ =γ∘1(z1)+γ∘2(z2).

This theorem allows us to easily compute the polars for and in terms of the polars of and , which are the dual norms, the spectral norm and infinity norm, respectively.

###### Corollary 4.2 (Explicit variational formulae for (($\text{SPCP}_{\text{sum}}$)) and (($\text{SPCP}_{\text{max}}$))).

We have

 ϕ∘sum(Z1,Z2) =max{|||Z1|||2,1λ% sum∥Z2∥∞} (4.10) ϕ∘max(Z1,Z2) Missing dimension or its units for \kern

where denotes the spectral norm (largest eigenvalue of ).

This result was also obtained by (van den Berg & Friedlander, 2011, Section 9), but is stated only for norms. Theorem 4.1 applies to gauges, and in particular now allows asymmetric gauges, so non-negativity constraints can be easily modeled.

We now have closed form solutions for in (4.8) for both formulations of interest. The remaining challenge is to design a fast solver for (4.3) for formulations (($\text{SPCP}_{\text{sum}}$)) and (($\text{SPCP}_{\text{max}}$)). We focus on this challenge in the remaining sections of the paper. We also discuss the advantage of (($\text{SPCP}_{\text{max}}$)) from this computational perspective.

## 5 Projections

In this section, we consider the computational issues of projecting onto the set defined by . For this is straightforward since the set is just the product set of the nuclear norm and norm balls, and efficient projectors onto these are known. In particular, projecting an matrix (without loss of generality let ) onto the nuclear norm ball takes operations, and projecting it onto the -ball can be done on operations using fast median-finding algorithms (Brucker, 1984; Duchi et al., 2008).

For , the projection is no longer straightforward. Nonetheless, the following lemma shows this projection can be efficiently implemented.

###### Proposition 5.1.

(van den Berg & Friedlander, 2011, Section 5.2) Projection onto the scaled -ball, that is, for some , can be done in time.

The proof of the proposition follows by noting that the solution can be written in a form depending only on a single scalar parameter, and this scalar can be found by sorting followed by appropriate summations. We conjecture that fast median-finding ideas could reduce this to in theory, the same as the optimal complexity for the -ball.

Armed with the above proposition, we state an important lemma below. For our purposes, we may think of as a vector in rather than a matrix in .

###### Lemma 5.2.

(van den Berg & Friedlander, 2011, Section 9.2) Let and , and let be any ordering of the elements of . Then the projection of onto the ball is , where is the projection onto the scaled -ball .

###### Sketch of proof.

We need to solve

 min{(L′,S′)∣ϕsum(L′,S′)≤1}12∥∥L′−L∥∥2F+12∥∥S′−S∥∥2F.

Alternatively, solve

 Missing dimension or its units for \kern

The inner minimization is equivalent to projecting onto the nuclear norm ball, and this is well-known to be soft-thresholding of the singular values. Since it depends only on the singular values, recombining the two minimization terms gives exactly a joint projection onto a scaled -ball. ∎

###### Remark 5.1.

All the references to the -ball can be replaced by the intersection of the -ball and the non-negative cone, and the projection is still efficient. As noted in Section 4, imposing non-negativity constraints is covered by the gauge results of Theorem 4.1 and Corollary 4.2. Therefore, both the variational and efficient computational framework can be applied to this interesting case.

## 6 Solving the Sub-Problem via Projected Quasi-Newton Methods

In order to accelerate the approach, we can use quasi-Newton (QN) methods since the objective has a simple structure.1 The main challenge here is that for the term, it is tricky to deal with a weighted quadratic term (whereas for , we can obtain a low-rank Hessian and solve it efficiently via coordinate descent).

We wish to solve (($\text{flip-SPCP}_{\text{max}}$)). Let be the full variable, so we can write the objective function as . To simplify the exposition, we take to be the matrix, but the presented approach applies to general linear operators (including terms like ). The matrix structure of and is not important here, so we can think of them as vectors instead of matrices.

The gradient is . For convenience, we use and

 ∇f(X)=(∇Lf(X)∇Sf(X))=AT(r(X)r(X)),rk≡r(Xk).

The Hessian is . We cannot simultaneously project onto their constraints with this Hessian scaling (doing so would solve the original problem!), since the Hessian removes separability. Instead, we use to approximate the cross-terms.

The true function is a quadratic, so the following quadratic expansion around is exact:

 f(L,S)=f(Xk) +⟨(∇Lf(Xk)∇Sf(Xk)),(L−LkS−Sk)⟩ +⟨(L−LkS−Sk),∇2f(L−LkS−Sk)⟩ =f(Xk) +⟨(rkrk),(L−LkS−Sk)⟩ +⟨(L−LkS−Sk),(1111)(L−LkS−Sk)⟩ =f(Xk) +⟨(rkrk),(L−LkS−Sk)⟩ +⟨(L−LkS−Sk),(L−Lk+S−SkL−Lk+S−Sk)⟩

The coupling of the second order terms, shown in bold, prevents direct 1-step minimization of , subject to the nuclear and 1-norm constraints. The FISTA (Beck & Teboulle, 2009) and spectral gradient methods (SPG) (Wright et al., 2009b) replace the Hessian with the upper bound , which solves the coupling issue, but potentially lose too much second order information. After comparing FISTA and SPG, we use the SPG method for solving (($\text{flip-SPCP}_{\text{sum}}$)). However, for (($\text{flip-SPCP}_{\text{max}}$)) (and for ((Lag-SPCP)), which has no constraints but rather non-smooth terms, which can be treated like constraints using proximity operators), the constraints are uncoupled and we can take a “middle road” approach, replacing

 ⟨(L−LkS−Sk),(L−Lk+S−SkL−Lk+S−Sk)⟩

with

 ⟨(L−LkS−Sk),(L−Lk+Sk−Sk−1Lk+1−Lk+S−Sk)⟩.

The first term is decoupled, allowing us to update , and then this is plugged into the second term in a Gauss-Seidel fashion. In practice, we also scale this second-order term with a number slightly greater than but less than (e.g., ) which leads to more robust behavior. We expect this “quasi-Newton” trick to do well when is similar to .

## 7 Numerical Results

The numerical experiments are done with the algorithms suggested in this paper as well as code from PSPG (Aybat et al., 2013), NSA (Aybat & Iyengar, 2014), and ASALM (Tao & Yuan, 2011)2. We modified the other software as needed for testing purposes. PSPG, NSA and ASALM all solve (($\text{SPCP}_{\text{sum}}$)), but ASALM has another variant which solves ((Lag-SPCP)) so we test this as well. All three programs also use versions of PROPACK from Larsen (1998); Becker & Candès (2008) to compute partial SVDs. Since the cost of a single iteration may vary among the solvers, we measure error as a function of time, not iterations. When a reference solution is available, we measure the (relative) error of a trial solution as . The benchmark is designed so the time required to calculate this error at each iteration does not factor into the reported times. Since picking stopping conditions is solver dependent, we show plots of error vs time, rather than list tables. All tests are done in Matlab and the dominant computational time was due to matrix multiplications for all algorithms; all code was run in the same quad-core  GHz i7 computer.

For our implementations of the (($\text{flip-SPCP}_{\text{max}}$)), (($\text{flip-SPCP}_{\text{sum}}$)) and ((Lag-SPCP)), we use a randomized SVD (Halko et al., 2011). Since the number of singular values needed is not known in advance, the partial SVD may be called several times (the same is true for PSPG, NSA and ASALM). Our code limits the number of singular values on the first two iterations in order to speed up calculation without affecting convergence. Unfortunately, the delicate projection involved in (($\text{flip-SPCP}_{\text{sum}}$)) makes incorporating a partial SVD to this setting more challenging, so we use Matlab’s dense SVD routine.

### 7.1 Synthetic test with exponential noise

We first provide a test with generated data. The observations with and were created by first sampling a rank matrix with random singular vectors (i.e., from the Haar measure) and singular values drawn from a uniform distribution with mean , and then adding exponential random noise (with mean equal to one tenth the median absolute value of the entries of ). This exponential noise, which has a longer tail than Gaussian noise, is expected to be captured partly by the term and partly by the term.

Given , the reference solution was generated by solving ((Lag-SPCP)) to very high accuracy; the values and were picked by hand tuning to find a value such that both and are non-zero. The advantage to solving ((Lag-SPCP)) is that knowledge of allows us to generate the parameters for all the other variants, and hence we can test different problem formulations.

With these parameters, was rank 17 with nuclear norm , had 54 non-zero entries (most of them positive) with norm , the normalized residual was , and , , , and .

Results are shown in Fig. 1. Our methods for (($\text{flip-SPCP}_{\text{max}}$)) and ((Lag-SPCP)) are extremely fast, because the simple nature of these formulations allows the quasi-Newton acceleration scheme of Section 6. In turn, since our method for solving (($\text{SPCP}_{\text{max}}$)) uses the variational framework of Section 4 to solve a sequence of (($\text{flip-SPCP}_{\text{max}}$)) problems, it is also competitive (shown in cyan in Figure 1). The jumps are due to re-starting the sub-problem solver with a new value of , generated according to (4.7).

Our proximal gradient method for (($\text{flip-SPCP}_{\text{sum}}$)), which makes use of the projection in Lemma 5.2, converges more slowly, since it is not easy to accelerate with the quasi-Newton scheme due to variable coupling, and it does not make use of fast SVDs. Our solver for (($\text{SPCP}_{\text{sum}}$)), which depends on a sequence of problems (($\text{flip-SPCP}_{\text{sum}}$)), converges slowly.

The ASALM performs reasonably well, which was unexpected since it was shown to be worse than NSA and PSPG in Aybat et al. (2013); Aybat & Iyengar (2014). The PSPG solver converges to the wrong answer, most likely due to a bad choice of the smoothing parameter ; we tried choosing several different values other than the default but did not see improvement for this test (for other tests, not shown, tweaking helped significantly). The NSA solver reaches moderate error quickly but stalls before finding a highly accurate solution.

### 7.2 Synthetic test from Aybat & Iyengar (2014)

We show some tests from the test setup of Aybat & Iyengar (2014) in the case. The default setting of was used, and then the NSA solver was run to high accuracy to obtain a reference solution . From the knowledge of , one can generate , but not and , and hence we did not test the solvers for ((Lag-SPCP)) in this experiment. The data was generated as , where was sampled by multiplication of and normal Gaussian matrices, had randomly chosen entries uniformly distributed within , and was white noise chosen to give a SNR of  dB. We show three tests that vary the rank from and the sparsity ranging from . Unlike Aybat & Iyengar (2014), who report error in terms of a true noiseless signal , we report the optimization error relative to .

For the first test (with and ), had rank and nuclear norm ; had of its elements nonzero and norm , and . The other parameters were , , , and . An interesting feature of this test is that while is low-rank, is nearly low-rank but with a small tail of significant singular values until number . We expect methods to converge quickly to low-accuracy where only a low-rank approximation is needed, and then slow down as they try to find a larger rank highly-accurate solution.

The results are shown in Fig. 2. Errors barely dip below (for comparison, an error of is achieved by setting ). The NSA and PSPG solvers do quite well. In contrast to the previous test, ASALM does poorly. Our methods for (($\text{flip-SPCP}_{\text{sum}}$)), and hence (($\text{SPCP}_{\text{sum}}$)), are not competitive, since they use dense SVDs. We imposed a time-limit of about one minute, so these methods only manage a single iteration or two. Our quasi-Newton method for (($\text{flip-SPCP}_{\text{max}}$)) does well initially, then takes a long time due to a long partial SVD computation. Interestingly, (($\text{SPCP}_{\text{max}}$)) does better than pure (($\text{flip-SPCP}_{\text{max}}$)). One possible explanation is that it chooses a fortuitous sequence of values, for which the corresponding (($\text{flip-SPCP}_{\text{max}}$)) subproblems become increasingly hard, and therefore benefit from the warm-start of the solution of the easier previous problem. This is consistent with empirical observations regarding continuation techniques, see e.g., (van den Berg & Friedlander, 2008; Wright et al., 2009b).

Figure 3 is the same test but with and , and the conclusions are largely similar.

### 7.3 Cloud removal

Figure 4 shows 15 images of size from the MODIS satellite,3 after some transformations to turn images from different spectral bands into one grayscale images. Each image is a photo of the same rural location but at different points in time over the course of a few months. The background changes slowly and the variability is due to changes in vegetation, snow cover, and different reflectance. There are also outlying sources of error, mainly due to clouds (e.g., major clouds in frames 5 and 7, smaller clouds in frames 9, 11 and 12), as well as artifacts of the CCD camera on the satellite (frame 4 and 6) and issues stitching together photos of the same scene (the lines in frames 8 and 10).

There are hundreds of applications for clean satellite imagery, so removing the outlying error is of great practical importance. Because of slow changing background and sparse errors, we can model the problem using the robust PCA approach. We use the (($\text{flip-SPCP}_{\text{max}}$)) version due to its speed, and pick parameters by using a Nelder-Mead simplex search. For an error metric to use in the parameter tuning, we remove frame from the data set (call it ) and set to be frames 2–15. From this training data , the algorithm generates and . Since is a matrix, it has far from full column span. Thus our error is the distance of from the span of , i.e., .

Our method takes about 11 iterations and 5 seconds, and uses a dense SVD instead of the randomized method due to the high aspect ratio of the matrix. Some results of the obtained outputs are in Fig. 5, where one can see that some of the anomalies in the original data frames are picked up by the term and removed from the term. Frame 4 has what appears to be a camera pixel error; frame 6 has another artificial error (that is, caused by the camera and not the scene); and frame 12 has cloud cover.

## 8 Conclusions

In this paper, we reviewed several formulations and algorithms for the RPCA problem. We introduced a new denoising formulation (($\text{SPCP}_{\text{max}}$)) to the ones previously considered, and discussed modeling and algorithmic advantages of denoising formulations (($\text{SPCP}_{\text{max}}$)) and (($\text{SPCP}_{\text{sum}}$)) compared to flipped versions (($\text{flip-SPCP}_{\text{max}}$)) and (($\text{flip-SPCP}_{\text{sum}}$)). In particular, we showed that these formulations can be linked using a variational framework, which can be exploited to solve denoising formulations using a sequence of flipped problems. For (($\text{flip-SPCP}_{\text{max}}$)), we proposed a quasi-Newton acceleration that is competitive with state of the art, and used this innovation to design a fast method for (($\text{SPCP}_{\text{max}}$)) through the variational framework. The new methods were compared against prior art on synthetic examples, and applied to a real world cloud removal application application using publicly available MODIS satellite data.

\pdfbookmark

[1]ReferencesrefSection

### Footnotes

1. We use “quasi-Newton” to mean an approximation to a Newton method and it should not be confused with methods like BFGS
2. PSPG, NSA and ASALM available from the experiment package at http://www2.ie.psu.edu/aybat/codes.html

### References

1. Aravkin, A. Y., Burke, J., and Friedlander, M. P. Variational properties of value functions. SIAM J. Optimization, 23(3):1689–1717, 2013a.
2. Aravkin, A. Y., Kumar, R., Mansour, H., Recht, B., and Herrmann, F. J. A robust SVD-free approach to matrix completion, with applications to interpolation of large scale data. 2013b.
3. Aybat, N., Goldfarb, D., and Ma, S. Efficient algorithms for robust and stable principal component pursuit. Computational Optimization and Applications, (accepted), 2013.
4. Aybat, N. S. and Iyengar, G. An alternating direction method with increasing penalty for stable principal component pursuit. Computational Optimization and Applications, (submitted), 2014.
5. Baldassarre, L., Cevher, V., McCoy, M., Tran Dinh, Q., and Asaei, A. Convexity in source separation: Models, geometry, and algorithms. Technical report, 2013.
6. Beck, A. and Teboulle, M. A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems. SIAM J. Imaging Sciences, 2(1):183–202, January 2009.
7. Becker, S. and Candès, E. Singular value thresholding toolbox, 2008. Available from http://svt.stanford.edu/.
8. Brucker, P. An O(n) algorithm for quadratic knapsack problems. Operations Res. Lett., 3(3):163 – 166, 1984.
9. Candès, E. J., Li, X., Ma, Y., and Wright, J. Robust principal component analysis? J. Assoc. Comput. Mach., 58(3):1–37, May 2011.
10. Chandrasekaran, V., Sanghavi, S., Parrilo, P. A., and Willsky, A. S. Sparse and low-rank matrix decompositions. In SYSID 2009, Saint-Malo, France, July 2009.
11. Chandrasekaran, V., Parrilo, P. A., and Willsky, A. S. Latent variable graphical model selection via convex optimization. Ann. Stat., 40(4):1935–2357, 2012.
12. Chen, Caihua, He, Bingsheng, Ye, Yinyu, and Yuan, Xiaoming. The direct extension of admm for multi-block convex minimization problems is not necessarily convergent. Optimization Online, 2013.
13. Combettes, P. L. and Pesquet, J.-C. A Douglas–Rachford Splitting Approach to Nonsmooth Convex Variational Signal Recovery. IEEE J. Sel. Topics Sig. Processing, 1(4):564–574, December 2007.
14. Duchi, J., Shalev-Shwartz, S., Singer, Y., and Chandra, T. Efficient projections onto the l1-ball for learning in high dimensions. In Intl. Conf. Machine Learning (ICML), pp. 272–279, New York, July 2008. ACM Press.
15. Ganesh, A., Lin, Z., Wright, J., Wu, L., Chen, M., and Ma, Y. Fast algorithms for recovering a corrupted low-rank matrix. In Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), pp. 213–215, Aruba, Dec. 2009.
16. Halko, N., Martinsson, P.-G., and Tropp, J. A. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM review, 53(2):217–288, 2011.
17. Huber, P. J. Robust Statistics. John Wiley and Sons, 2 edition, 2004.
18. Kyrillidis, Anastasios and Cevher, Volkan. Matrix recipes for hard thresholding methods. Journal of Mathematical Imaging and Vision, 48(2):235–265, 2014.
19. Larsen, R. M. Lanczos bidiagonalization with partial reorthogonalization. Tech. Report. DAIMI PB-357, Department of Computer Science, Aarhus University, September 1998.
20. Lee, J., Recht, B., Salakhutdinov, R., Srebro, N., and Tropp, J.A. Practical large-scale optimization for max-norm regularization. In Neural Information Processing Systems (NIPS), Vancouver, 2010.
21. Lin, Z., Chen, M., and Ma, Y. The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices. arXiv preprint arXiv:1009.5055, 2010.
22. Mansour, H. and Vetro, A. Video background subtraction using semi-supervised robust matrix completion. In To appear in IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 2014.
23. Peng, Y., Ganesh, A., Wright, J., Xu, W., and Ma, Y. RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated images. IEEE Trans. Pattern Analysis and Machine Intelligence, 34(11):2233–2246, 2012.
24. Recht, Benjamin and Ré, Christopher. Parallel stochastic gradient algorithms for large-scale matrix completion. Math. Prog. Comput., pp. 1–26, 2011.
25. Rockafellar, R. T. Convex analysis. Princeton Mathematical Series, No. 28. Princeton University Press, Princeton, N.J., 1970.
26. Shen, Y., Wen, Z., and Zhang, Y. Augmented Lagrangian alternating direction method for matrix separation based on low-rank factorization. Optimization Methods and Software, 29(2):239–263, March 2014.
27. Tao, M. and Yuan, X. Recovering low-rank and sparse components of matrices from incomplete and noisy observations. SIAM J. Optimization, 21:57–81, 2011.
28. van den Berg, E. and Friedlander, M. P. Probing the Pareto frontier for basis pursuit solutions. SIAM J. Sci. Computing, 31(2):890–912, 2008. software: http://www.cs.ubc.ca/~mpf/spgl1/.
29. van den Berg, E. and Friedlander, M. P. Sparse optimization with least-squares constraints. SIAM J. Optimization, 21(4):1201–1229, 2011.
30. Wright, J., Ganesh, A., Rao, S., and Ma, Y. Robust principal component analysis: Exact recovery of corrupted low-rank matrices by convex optimization. In Neural Information Processing Systems (NIPS), 2009a.
31. Wright, S. J., Nowak, R. D., and Figueiredo, M. A. T. Sparse Reconstruction by Separable Approximation. IEEE Trans. Sig. Processing, 57(7):2479–2493, July 2009b.
32. Zhang, Z., Liang, X., Ganesh, A., and Ma, Y. TILT: Transform invariant low-rank textures. In Kimmel, R., Klette, R., and Sugimoto, A. (eds.), Computer Vision – ACCV 2010, volume 6494 of Lecture Notes in Computer Science, pp. 314–328. Springer, 2011.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters


text:
# The Pantheon+ Analysis: Forward-Modeling the Dust and Intrinsic Colour Distributions of Type Ia Supernovae, and Quantifying their Impact on Cosmological Inferences

Recent studies have shown that the observed colour distributions of Type Ia SNe (SNIa) are well-described by a combination of distributions from dust and intrinsic colour. Here we present a new forward-modeling fitting method (Dust2Dust) to measure the parent dust and colour distributions, including their dependence on host-galaxy mass. At each fit step, the SNIa selection efficiency is determined from a large simulated sample that is re-weighted to reflect the proposed distributions. We use five separate metrics to constrain the Dust2Dust parameters: distribution of fitted light-curve colour $c$, cosmological residual trends with $c$, cosmological residual scatter with $c$, fitted colour-luminosity relationship $\beta_{\rm SALT2}$, and intrinsic scatter $\sigma_{\rm int}$. Using the Pantheon+ data sample, we present results for a Dust2Dust fit that includes 4 parameters describing intrinsic colour variations and 8 parameters describing dust. Furthermore, we propagate the Dust2Dust parameter uncertainties and covariance to the dark energy equation-of-state $w$ and Hubble constant H$_0$: we find $\sigma_w = 0.005$ and $\sigma_{\textrm{H}_0} = 0.145~$km/s/Mpc. The Dust2Dust code is publically available.

### Cited Authors

• Popovic, B; Brout, D; Kessler, R; Scolnic, D

### Published Date

• December 10, 2021


text:
# mathlibdocumentation

algebraic_geometry.Scheme

# The category of schemes #

A scheme is a locally ringed space such that every point is contained in some open set where there is an isomorphism of presheaves between the restriction to that open set, and the structure sheaf of Spec R, for some commutative ring R.

A morphism of schemes is just a morphism of the underlying locally ringed spaces.

structure algebraic_geometry.Scheme  :
Type (max (u_1+1) (u_2+1))

We define Scheme as a X : LocallyRingedSpace, along with a proof that every point has an open neighbourhood U so that that the restriction of X to U is isomorphic, as a space with a presheaf of commutative rings, to Spec.PresheafedSpace R for some R : CommRing.

(Note we're not asking in the definition that this is an isomorphism as locally ringed spaces, although that is a consequence.)

Every Scheme is a LocallyRingedSpace.

Equations
@[instance]

Schemes are a full subcategory of locally ringed spaces.

Equations

The spectrum of a commutative ring, as a scheme.

Equations

The induced map of a ring homomorphism on the ring spectra, as a morphism of schemes.

Equations
@[simp]
theorem algebraic_geometry.Scheme.Spec_map_comp {R S T : CommRing} (f : R S) (g : S T) :
@[simp]

The spectrum, as a contravariant functor from commutative rings to schemes.

Equations
@[simp]

The empty scheme, as Spec 0.

Equations
@[instance]
Equations
@[instance]
Equations

The global sections, notated Gamma.

Equations
@[simp]


text:
# How do I find the area inside the cardioid r=1+costheta?

Aug 8, 2018

$4 \pi$

#### Explanation:

We will integrate the area differential for a polar function:
$\mathrm{dA} = d \left(\frac{1}{2} {r}^{2} \theta\right) = r \theta d r + \frac{1}{2} {r}^{2} d \theta$

We can integrate this just $d \theta$:
$\mathrm{dA} = \left(r \theta \frac{\mathrm{dr}}{d \theta} + \frac{1}{2} {r}^{2}\right) d \theta$

This yields
$A = \int \setminus \setminus \mathrm{dA} = {\int}_{0}^{2 \pi} \left(\left[1 + \cos \theta\right] \cdot \theta \cdot \left(- \sin \theta\right) + \frac{1}{2} {\left(1 + \cos \theta\right)}^{2}\right) d \theta$
$A = {\int}_{0}^{2 \pi} \left[\frac{3}{4} + \frac{1}{4} \cos 2 \theta + \cos \theta - \theta \sin \theta - \frac{1}{2} \theta \sin 2 \theta\right] d \theta$

This integral is pretty simple, with a few integrations by parts.

The ultimate answer turns up to be $A = 4 \pi$.