text:
# On 2D Inverse Problems/Kernel of Dirichlet-to-Neumann map

The continuous analog of the matrix representation of a Dirichlet-to-Neumann operator for a domain $\Omega$ is its kernel. It is a distribution defined on the Cartesian product of the boundary of the domain w/itself, such that if

$\Lambda f = g,$

then

$g(\phi) = \int_{\partial\Omega}K(\phi,\theta)f(\theta)d\theta,$

where $\phi$ and $\theta$ parametrise the boundary.

For the case of a half-plane in two-dimensions with uniform conductivity 1 one can calculate the kernel explicitly. Because the domain in consideration is shift invariant, the kernel is a convolution.

$K(\phi,\theta) = k(\phi-\theta),$

where k is a distribution on a line. Therefore, the calculation reduces to solving the Dirichlet problem for a $\delta_0$-function at the origin and taking normal derivative at the boundary line.

Exercise (**). Complete the calculation of the kernel K for the half-plane to show that:

$K(x,y) = \frac{-1}{\pi(x-y)^2}$

off the diagonal.

Exercise (*). Prove that for rotation invariant domain (disc w/ conductivity depending only on radius) the kernel of Dirichlet-to-Neumann map is a convolution.

The Hilbert transform gives a correspondence between boundary values of harmonic function and its harmonic conjugate.

$H:u|_{\partial\Omega}\rightarrow v|_{\partial\Omega},$

where

$f(z) = u(z) + iv(z)$

is an analytic function in the domain.

Exercise (*). Prove that for the case of the complex half-plane C+ the Hilbert transform is given by the following formula:

$H_{C^+}f(y) =\frac{1}{\pi} \ \text{p.v.} \int_{-\infty}^{\infty} \frac{f(x)}{y-x}dx.$

Exercise (*). Differentiate under integral sign the formula above to obtain the kernel representation for the Dirichlet-to-Neumann operator for the uniform half plane.

To define discrete Hilbert transform for a planar network, one considers the corresponding graph together w/its dual.


text:
# Difference between revisions of "1994 AIME Problems/Problem 13"

## Problem

The equation

$x^{10}+(13x-1)^{10}=0\,$

has 10 complex roots $r_1, \overline{r_1}, r_2, \overline{r_2}, r_3, \overline{r_3}, r_4, \overline{r_4}, r_5, \overline{r_5},\,$ where the bar denotes complex conjugation. Find the value of

$\frac 1{r_1\overline{r_1}}+\frac 1{r_2\overline{r_2}}+\frac 1{r_3\overline{r_3}}+\frac 1{r_4\overline{r_4}}+\frac 1{r_5\overline{r_5}}.$

format: plain


text:
# Complementary learning for AI-based predictive quality and maintenance

Using machine learning, deep learning, and cognitive computing in concert can help enterprises gain competitive edges.

May 3, 2018
Fractal kaleidoscope (source: Pixabay)

The cost of unplanned downtime due to equipment failure across all industries is huge. Back in 2014, Aberdeen estimated the cost to businesses to be $164,000 an hour, on average. By 2016, that figure had skyrocketed by 60% to$260,000 an hour. Monetary losses aside, the impacts of unplanned outages for a business can be devastating—according to a 2017 study by Vanson-Bourne, 46% of businesses said they were not able to deliver services to customers and 37% said they lost production time on a critical piece of equipment or other asset.

Take the airline industry as an example. Flight delays due to mechanical problems are common. Sometimes the issue gets fixed in timely manner. But sometimes it takes hours, and, in worst-case scenarios, flights are canceled. The costs to the passengers can be high—missed connections, missed meetings, missed vacations, and more. The cost to the airline is immense: lost revenue, lost customer satisfaction, and lost brand value in addition to all the money—labor and parts—it takes to fix the problem. Consider the overtime, the shipping costs of rushed parts, and all the lost revenue, and you see how it adds up. In fact, unplanned maintenance of equipment costs three to nine times more than planned maintenance.

## Learn faster. Dig deeper. See farther.

Join the O'Reilly online learning platform. Get a free trial today and find answers on the fly, or master something new and useful.

Clearly, a lot is at stake. Which is why many industries that focus on product quality and maintenance are turning to next-generation artificial intelligence (AI)-based predictive quality and maintenance (PQM) solutions.

In this article, I’ll go over what PQM solutions are and explain how they are being made dramatically more effective by leveraging AI. I’ll also explore complementary learning and look at how today’s leading AI-based PQM solutions incorporate it. Finally, I’ll give a real-world example of AI-based PQM in action.

## What are PQM solutions?

PQM solutions, which harness data gathered by both the Internet of Things (IoT) and data from traditional legacy systems, focus on detecting and addressing quality and maintenance issues before they turn into serious problems—for example, problems that can cause unplanned downtime.

Unplanned downtime is a major cost driver in any industry that must maintain large inventories of capital assets. For an airline, for example, delaying flights due to unplanned maintenance can cost thousands of dollars each minute. Unplanned shutdowns of oil platforms can run into the millions of dollars. And in manufacturing plants, the costs of disruptions go directly to the bottom line. It is the goal of every organization to eliminate unplanned downtime in favor of planned maintenance. PQM solutions can help with planned maintenance, too, by shortening maintenance operations windows.

PQM solutions help businesses prioritize how to best allocate scarce resources, resolve issues, and plan ahead, thereby keeping—or achieving—a competitive edge. Recent research suggests that the market for PQM solutions will grow from $2.2 billion in 2017 to$10.9 billion by 2022, a 39% annual growth rate. And, of the top 10 uses predicted for AI in 2025, PQM comes in fourth place.

## Artificial intelligence enters the picture

Traditionally, PQM solutions crunched numbers and came up with average statistics to predict when quality corrections or maintenance were required.

Now, with the availability of much bigger data sets and the latest developments in AI, it’s a whole new game. AI-based PQM solutions differ from traditional predictive quality and maintenance solutions because they analyze the actual condition of a product rather than just using average or expected statistics to predict when quality corrections or maintenance might be required.

AI-based PQM solutions use several technologies in concert, including machine learning, deep learning, and cognitive computing:

• Machine learning: focuses on real-world problems by processing—and learning from—large amounts of data
• Deep learning: uses neural networks to be able to sort through nearly unimaginable volumes of data to come to conclusions
• Cognitive computing: a subset of AI that attempts to mimic the way humans think. And a very important subset of cognitive computing is associative-memory learning and reasoning, which mimics the way humans learn, remember, and reason by making associations

## Complementary learning: The future of PQM solutions

Because each type of AI is good at solving different problems, applying them simultaneously is the key to success. Complementary learning in the context of PQM applications involves combining all these types of AI—machine learning, deep learning, and cognitive computing—to get insight into quality and maintenance issues.

In effect, a PQM solution that embraces complementary learning first uses machine learning and deep learning to answer the question, “What is the problem?” Then, cognitive computing answers such questions as: “Have I ever seen this before? What type of a problem is it? Who knows how to fix this? What caused this problem? And will it happen again?”

## AI-based PQM solution in action

Accenture is a global professional services company that provides a broad range of services and solutions in business strategy, consulting, digital, technology, and operations. One of its many services focuses on software testing.

This is a rich area for AI-based PQM solutions for a couple of reasons. First, many organizations have not changed their software-testing techniques in the past 20 years. Thus, software testing often involves a significant amount of manual labor. For some mission-critical systems, test engineers spend up to 90% of their time managing test cases and documenting them rather than actually testing.

Second, the overall quality of software has traditionally been treated as just one isolated component in the software development life cycle. For example, QA used to be considered a step outside of the product development cycle—involve QA when everything is done. But testing can no longer be performed in silos outside of the development process. Now, with devices and connections getting more complicated, QA needs to be integrated as part of the software’s design and production process. If QA sees a trend that a certain part or certain function causes frequent issues, they can report such a trend to the product team to find the root cause and fix.

Accenture’s platform uses Intel Saffron AI in its software testing solution to accelerate automation, spot trends, manage risks, and continuously respond to customer feedback. According to Accenture, this approach has reduced time-to-market and cost of testing by more than 20%, and can also save 30% to 50% of time spent over-engineering.

“Testing is transforming into quality engineering, where artificial intelligence-driven analytics is at the core of driving productivity and agility,” says Kishore Durg, senior managing director and global testing lead for Accenture. “That is what the Accenture Touchless Testing Platform is designed to do.”


text:
# Robotics

a syllabus

## Course description

The course begins with a survey of the evolving field of robotics, including major categories such as mobile and manipulation robotics. Topics surveyed include effectors, actuators, locomotion, manipulation, sensors, feedback control, control architectures, navigation, emergence, and group robotics. Upon conclusion of the survey, the focus of the course turns to a special topic chosen by the instructor. The instructor could choose, for instance, swarm robotics or robot kinematics and dynamics—or a project. (Adopted from course catalog.)

The current course special topic is robot planning and kinematics. We will use the page
Minicourse on Robot Planning and Kinematics
as a secondary syllabus.

## General information

Actual office hours (CH 103C)
T 4-5, W 1-2, Th 4-5
Virtual office hours (zoom link, make appointment!)
F 1-3
Virtual office hours appointments
make appointment
Time (synchronous)
W 5-6:50
Location
Zoom!
https://us02web.zoom.us/j/89025317510
unlisted playlist
Website
ricopic.one/courses/robotics
Minicourse website
ricopic.one/courses/robotics-planning-kinematics
Moodle
Robotics Moodle

## Textbooks

For textbooks, see the Minicourse on Robot Planning and Kinematics.

## Robotics bük

A digital text with integrated video lectures (and more), called a bük, I’m writing is available robotics.ricopic.one (henceforth: Ro).

To get familiar with the use of a bük, see this bük.

## Videos

Asynchronous video lectures are integrated into Ro and any related minicourse material.

A synchronous session will be held each week on Zoom (see General information).

Synchronous sessions are recorded and will be available within about 24 hours at the unlisted YouTube playlist shown in General information.

## Microsoft Teams

Everyone is required to join the Microsoft Teams team MME 565. We’ll use it to communicate with each other during the semester. Join here MME 565. That’s a signup link. Be sure to join the channel General.

## Schedule

The following schedule is tentative. The rest of the schedule can be found at the page Minicourse on Robot Planning and Kinematics.

day topics week bük due
Introduction 1 bük bük
Ro [in]
Robot mechanics and control architectures 2 Ro [mec]
Ro [co]
See Minicourse on Planning and Kinematics 3+

## Assignments

Assignments are due late Saturday of the week they appear on the schedule. They should be turned in on Moodle.

For more, see the page Minicourse on Robot Planning and Kinematics.

Total grades in the course may be curved. They will be available on Moodle throughout the semester.

Assignments
100%

A grading rubric visible in Moodle will be used to grade assignments.

## Correlation of course & program outcomes

In keeping with the standards of the Department of Mechanical Engineering, each course is evaluated in terms of its desired outcomes and how these support the desired program outcomes. The following sections document the evaluation of this course.

### Desired course outcomes

Upon completion of the course, the following course outcomes are desired:

1. students will demonstrate the ability to use the fundamentals of robotics;
2. students will demonstrate the ability to analyze in detail systems of at least one contemporary topic in robotics;
3. students will demonstrate the ability to analyze robot actuators;
4. students will demonstrate the ability to analyze and design robot navigation; and
5. students will demonstrate the ability to analyze emergence.

### Desired program outcomes

In accordance with ABET’s student outcomes, our desired program outcomes are that mechanical engineering graduates have:

1. an ability to identify, formulate, and solve complex engineering problems by applying principles of engineering, science, and mathematics
2. an ability to apply engineering design to produce solutions that meet specified needs with consideration of public health, safety, and welfare, as well as global, cultural, social, environmental, and economic factors
3. an ability to communicate effectively with a range of audiences
4. an ability to recognize ethical and professional responsibilities in engineering situations and make informed judgments, which must consider the impact of engineering solutions in global, economic, environmental, and societal contexts
5. an ability to function effectively on a team whose members together provide leadership, create a collaborative and inclusive environment, establish goals, plan tasks, and meet objectives
6. an ability to develop and conduct appropriate experimentation, analyze and interpret data, and use engineering judgment to draw conclusions
7. an ability to acquire and apply new knowledge as needed, using appropriate learning strategies.

### Correlation of outcomes

The following table correlates the desired course outcomes with the desired program outcomes they support.

desired program outcomes
1 2 3 4 5 6 7
desired course outcomes 1

Cheating or plagiarism of any kind is not tolerated and will result in a failing grade (“F”) in the course. I take this very seriously. Engineering is an academic and professional discipline that requires integrity. I expect students to consider their integrity of conduct to be their highest consideration with regard to the course material.

Cheating is academic dishonesty as well as unprofessional for prospective teachers. Do not copy other students' assignments, have someone else write your papers or plagiarize published or unpublished materials, or submit work previously graded by other instructors. See Saint Martin's University Student Handbook. Students will be graded not only on their academic success, but on professional conduct as well. Students who fail to show professionalism in their academic or personal conduct (e.g. constant tardiness, excessive absences, and/or other unprofessional behavior) may earn a lower letter grade than the total of semester accumulated points, or may even earn a failing grade.

What is academic integrity? Saint Martin’s University is a community of faculty, students and staff engaged in the exchange of ideas in the ongoing pursuit of academic excellence. Essential to our mission is a focused commitment to scholarly values and intellectual integrity, and a respect for the ideas, beliefs and work of others. This commitment extends to all aspects of academic performance. All members are expected to abide by ethical standards both in their conduct and their exercise of responsibility to themselves and toward other members of the community. As an expression of our shared belief in the Benedictine tradition, we support the intellectual, social, emotional, physical and spiritual nurturing of students.

What is academic dishonesty? Saint Martin’s University defines academic dishonesty as violating the academic integrity of an assignment, test and/or evaluation of any coursework. This dishonest practice occurs when you seek to gain for yourself or another an academic advantage by deception or other dishonest means. You have a responsibility to understand the requirements that apply to particular assessments and to be aware of acceptable academic practice regarding the use of material prepared by others. Therefore, it is your responsibility to be familiar with the policies surrounding academic dishonesty as these may differ from other institutions.

## Access and accommodations

Your experience in this class is important to me. If you have already established accommodations with Disability Support Services (DSS), please communicate your approved accommodations to me at your earliest convenience so we can discuss your needs in this course.

If you have not yet established services through DSS, but have a temporary health condition or permanent disability that requires accommodations (conditions include but are not limited to mental health, attention-related, learning, vision, hearing, physical or health impacts), you are welcome to contact DSS at 360-438-4580 or dss.testing@stmartin.edu or smu.dss@stmartin.edu. DSS offers resources and coordinates reasonable accommodations for students with disabilities and/or temporary health conditions. Reasonable accommodations are established through an interactive process between you, your instructor(s) and DSS. It is the policy and practice of the Saint Martin’s University to create inclusive and accessible learning environments consistent with federal and state law.

## Sexual misconduct/sexual harassment reporting

Saint Martin’s University is committed to providing an environment free from sex discrimination, including sexual harassment and sexual violence. There are Title IX/sexual harassment posters around campus that include the contact information for confidential reporting and formal reporting. Confidential reporting is where you can talk about incidents of sexual harassment and gender-based crimes including sexual assault, stalking, and domestic/relationship violence. This confidential resource can help you without having to report your situation to the formal reporting process through the Dean of Students – Ms. Melanie Richardson, Associate VP of Human Resources – Ms. Cynthia Johnson, Public Safety – Ms. Sharon Schnebly, or the Office of the Provost – Dr. Kathleen Boyle, unless you request that they make a report. Please be aware that, in compliance with Title IX and under the Saint Martin’s University policies, educators must report incidents of sexual harassment and gender-based crimes including sexual assault, stalking, and domestic/relationship violence. If you disclose any of these situations in class, on papers, or to me personally, I am required to report it.

## University sanctioned activities

If a student is absent from class due to university sanctioned activities, such as sports, it is the student's responsibility to request that the absence be excused, otherwise, the absence will be recorded as unexcused. Absent students are responsible for catching up with the class, and if any assignments are due on the day of the absence, it is the student's responsibility to turn in the assignments on time (prior to class).

## Religious Accommodation

Saint Martin’s University, in honor of the sacredness of the individual, and being deeply rooted in the Catholic Benedictine tradition of higher education, values the many religious and spiritual practices of our campus community. Saint Martin’s University supports our students in their ongoing journey of becoming. In compliance with Washington State Law RCW 28B.137.010, Saint Martin’s University reasonably accommodates students for reasons of religious observances.


text:
# Machine mathematical reasoning by clever substitutions, How to do with AI

I have three equations that relates five variables {a, b, c, r, s} with a sum and two ratios.

Eq. 1: a = b + c;
Eq. 2: s = b / a;
Eq. 3: r = b / c.

Given two values for any of the five variables I get a solution. But, this is not the automation problem I want to solve.

I can have the solution of variable r by simply knowing s. This is solved by a "human algorithm" as follows.

1. Substitute a of Eq. 1 in Eq. 2.
2. Divide the second term of the new Eq. 2 by the variable c.
3. Replace b/c by the expression of Eq. 3.

That means s = r / (r+1).

The questions is -How can an AI algorithm solve this?, i.e. the machine should recognize that given the variable r she can obtain directly the variable s and do no require another variable.

What you're asking here is actually not reasoning but rather symbolic computation. Sympy is an example of a library which does it.

Explaining computer algebra systems could be a series of lectures and this is too broad for a single question. To give you an idea: (1) represent equations symbolically (2) bring them in a normal form (3) substitute.

Although it is not elegant, you can simplify the problem by brute Force: try all combinations of solving for a variable and substituting it in the other equations.

Related:

I would turn this problem into graph search in which vertices are "equivalent sets of equations" and edges are "valid operations" such as such as substitution, division and so on, which always ensure that all states are mathematically equivalent. A very simple example of a few states and edges (to requote from your question):

State 0:
Eq. 1: a = b + c;
Eq. 2: s = b / a;
Eq. 3: r = b / c.

Operation 0-->1:
Multiply Eq. 3 by c

State 1:
Eq. 1: a = b + c;
Eq. 2: s = b / a;
Eq. 3: rc = b.

Operation 0-->2:
Substitute Eq 1 into 2

State 2:
Eq. 1: a = b + c;
Eq. 2: s = b / (b+c);
Eq. 3: r = b / c.

Operation 2-->3:
Rearrange Eq 3 so B is on right hand side.

etc etc

Depending on how many operations you have, this does mean that the graph could be quite large, so you may not want to generate the entire graph (actually, the graph will be infinite if the number of operations allowed is unlimited). It's also possible that the graph will contain cycles, since there may be multiple ways to reach the same state.

You will have to think about how to develop data structures for representing the states and operations. Since this is a search problem, it will help to define what your "goal" states look like (e.g. states where s is defined in terms of r).

To search the graph, I recommend one of these approaches:

• Breadth-first search: this will find the shortest number of operations needed to reach any goal states
• Best-first search: if the graph is too large for breadth first, you may need to guide the search to the goal states using a heuristic, e.g. you could rate states by how few variables are needed on the LHS of equations or something like that -- some experimentation will be needed here.

The exact choice of search algorithm will depend on how big the graph is, which in turns depends on how many operations/edges you want.


text:
# How do I add a copyToClipboard function in the following code?

I have tried many ways to do it please help me with this. There’s a icon which gives a clipboard of the link to email. I have to remove that clipboard and add copytoclipboard when clicked on the icon. I have commented the previous code where it gave the clipboard of the link.

$scope.copyToClipBoard =function (tranxNo,messageID,userCode,docCode,processCode, origin,company,originCompany,isInterInstance) { // need messageID,us_login,InstanceId,Flg to mark this unique url. var userCode=$('.loggedInUserCode').val();
var curTrnxNo = $("#trnxNo").val(); var curBpCode =$("#bpCode").val();

var clipBoardId="";
if(docCode.indexOf("PR0")==0)
clipBoardId = messageID+"_"+curBpCode;
else
clipBoardId = curTrnxNo+"_"+curBpCode;
var clipBoardDiv=$(".clipBoardDiv_"+clipBoardId);$(".clipBoardInput_"+clipBoardId).val("");

var parmUrl="?callback=JSON_CALLBACK&authToken="+$scope.authToken +"&messageID="+messageID +"&processCode="+processCode +"&tranxNo="+tranxNo +"&userCode="+userCode +"&origin="+origin +"&company="+company +"&originCompany="+originCompany +"&isInterInstance="+isInterInstance; var mailURL=$scope.realPath+streamMailPath+parmUrl;
\$http.jsonp(mailURL).success(function(data) {
if(data!="undefined")
{
replyTo=data.reply_To;

//openClipBoardDialog(tranxNo,messageID,docCode,processCode,replyTo);

}
}).error(function(data, status, headers, config) {
});
};

Source: AngularJS


text:
1887

# News Picks

Physics Today’s online staff summarize the most important and interesting news about science from the world's top media outlets.

There are 70 posts for the selected month (July 2015).

### Filter by date

July 27, 2015 2:13 PM

### Massive black hole found during study of archival data

Scientific Computing: A black hole 100 million solar masses in size has been spotted in the process of devouring a star. While examining spectra of objects taken by the Sloan Digital Sky Survey, Andrea Merloni of the Max Planck Institute for Extraterrestrial Physics in Garching, Germany, and his colleagues found an unusual change between 1998 and 2005 in the spectrum of one galaxy located some 3.5 billion light-years away. Further confirmation of the high-energy emission came when the team discovered that both the European Space Agency's XMM-Newton and NASA's Chandra X-Ray Observatory, which had taken pictures of the region at about the same time, had also seen the massive flare. An analysis of the data from the three archival sources confirmed that the flare fit the model of a star being pulled into a massive black hole.

July 27, 2015 1:03 PM

### First data from new LHC run are presented in Vienna

Guardian: Researchers from CERN's Large Hadron Collider (LHC) are presenting their latest results at the European Physical Society's Conference on High Energy Physics in Vienna. The ATLAS team confirmed an unexpected finding from the CMS team during the previous LHC run—that regardless of the incident collision angle of the particles, there is an increased probability they are being emitted at a similar azimuthal angle. That behavior is expected in heavy particle collisions but not in proton–proton collisions. From the new run, ATLAS also released the first measurement of the jets of hadrons produced from the collision of quarks and gluons. The measurement, which is in line with the measurements taken during the previous run, reveals that the higher collision energy is producing noticeably more frequent jets.

July 27, 2015 11:45 AM

### Brain revealed through photoacoustic imaging

NPR: A new technique that combines “the speed and precision of light with the penetrating ability of sound” has been used to create high-speed, detailed three-dimensional images of a living mouse brain, writes Jon Hamilton for NPR. Developed by Lihong Wang of Washington University in St Louis, photoacoustic imaging uses a laser to send light pulses through the mouse’s skull and into its brain, where they bounce around and cause the molecules to vibrate. The vibrations emit distinctive sound waves, which can be used to monitor brain activity and study individual brain cells. The technique could one day be applied to other areas of the body to look for anomalies such as tumors in the breast, skin, and even individual blood cells.
July 27, 2015 11:40 AM

### "Bubble" greenhouse proposed for remote, arid regions

SciDev.Net: A sealed greenhouse that uses evaporation and condensation to desalinate saltwater for irrigation is being developed by researchers in Australia. In the design the evaporation equipment is located outside the greenhouse structure. Air is bubbled through water-filled columns to increase the air–water interface and encourage evaporation. Although the novel humidification–dehumidification process requires more energy than other desalination methods, such as solar-powered reverse osmosis, its simplicity and small scale could make it useful in isolated areas where local people with little training could set it up and run it. In addition, such a sealed structure would provide a cool, humid environment, which promotes plant growth, and would protect crops from insects and disease.
July 24, 2015 12:10 PM

### Mammoth DNA suggests climate change was primary driver in extinction

Science: Mammoths and many other large mammals suffered a massive population decline between 30 000 and 20 000 years ago that eventually led to their extinction. Human expansion, and the hunting associated with it, is widely considered to be a major factor in that decline, but a new study of DNA collected from mammoth fossils suggests that the populations were shrinking well before human expansion was significant. Alan Cooper of the University of Adelaide in Australia and his colleagues established a connection between the diversity of DNA found at specific sites and the size of the species population. Creating a timeline of temperature changes from sediment cores, they found that drops in DNA diversity, and therefore overall population size, were tied to periods of significant short-term warming.

July 24, 2015 11:34 AM

### Chiral silicon could be used in light-based computing

Science: Some schemes for photonic computing use light's polarization to store and transfer data. But if the data are to be transferred to a regular computer, a silicon-based means of measuring spin is needed, which so far has proved elusive to realize. Now Ritesh Agarwal of the University of Pennsylvania and his colleagues have discovered that by altering the structure of a silicon crystal, the material is no longer indistinguishable from its mirror image. That property allows silicon to be used to detect both polarization and electron spin because of the way light and electrons interact with the crystal. The availability of silicon detectors will significantly increase the rate at which photonic computing can transmit data.

July 23, 2015 3:15 PM

### Berkeley passes controversial cell-phone radiation ordinance

New York Times: The city of Berkeley, California, has passed a measure requiring cell-phone retailers to warn customers that prolonged use of the devices may be hazardous to their health. According to the Right to Know ordinance, people who regularly carry cell phones in a pocket or bra “may exceed the federal guidelines for exposure” to radio frequency radiation. However, the claim has not been substantiated by any scientific study. Furthermore, no adverse effects associated with low levels of radiation have been seen over the past 50 years and there has not been any noticeable increase in brain cancers, according to Jerrold Bushberg, a medical physicist and a professor of radiology and radiation oncology at the University of California, Davis, and a representative of the American Association of Physicists in Medicine. Regardless, proponents maintain that if there’s any possibility of risk, no matter how slight, the public should be made aware. The cell-phone industry has challenged the ordinance, and a federal hearing has been set for 6 August.
July 23, 2015 2:12 PM

Nature: In December, Tim Peake will become the first astronaut to visit the International Space Station (ISS) as part of the UK's space program. The milestone will mark a major step for the UK Space Agency, which was established in 2010, and exemplifies the country's interest in participating in human spaceflight. A long-time contributor to the European Space Agency (ESA), the UK has participated primarily on telescope and robotic probe projects. However, in 2012 the country began expanding its involvement in human spaceflight by investing €20 million ($22 million) in the ISS and another €16 million in the ESA's European Programme for Life and Physical Sciences, which studies the effects of spaceflight on humans onboard the ISS. July 23, 2015 12:40 PM ### Slips in ITER schedule mean China's fusion project may finish first New Scientist: ITER, the international prototype fusion energy reactor being built in France, has seen costs triple to$20 billion and has experienced construction delays that have pushed back completion to 2019. Designed to be the world's largest fusion facility, ITER is expected to be the first to produce more energy than was put into it. However, the delays may prevent it from achieving self-sustaining fusion reactions until the 2030s. That means that China, which is also participating in ITER, will likely finish construction on its own fusion project in time to reach that threshold first. The China Fusion Engineering Test Reactor is still in the design stage, but the country is heavily investing in the project and is likely learning quite a lot from its participation in ITER.

July 23, 2015 11:55 AM

### Foldable bridge created for use in emergencies

BBC: A new bridge design is being developed for use during natural disasters and other emergencies when local infrastructure has been damaged or is otherwise unusable. Called Mobile Bridge version 4.0, the scissors-like device can fold up compactly for transport and be readily expanded as needed. It’s the brainchild of Ichiro Ario of Hiroshima University and his colleagues, who recently tested a prototype at the Hongo River in Fukuyama, Japan. It can be fully deployed in just one hour and requires only a few people to set it up. The researchers say it is “the world’s fastest, largest, strongest, and lightest expanding temporary bridge.”
July 22, 2015 2:40 PM

### Lizard’s skin channels water to mouth

BBC: Certain lizards living in arid regions have been shown to collect water via a sophisticated capillary system running between their scales. Researchers studying the Texas horned lizard discovered that the interlinked capillary channels, which are partially enclosed by the lizard's overlapping scales, are wider in the direction of the tail and grow narrower in the direction of the mouth. Water entering the system is squeezed forward because of the channels' shape and width. Now Philipp Comanns of Aachen University and his colleagues have created a laboratory prototype of the "passive, directional liquid transport" by laser etching channels into a glass-like plastic and leaving raised “scales” in between. Such a capillary system could have practical applications in distilleries, heat exchangers, and small medical devices.
July 22, 2015 11:40 AM

### Explosion at NIST similar to a meth lab blast

Washington Post: Over the weekend, an explosion at NIST in Gaithersburg, Maryland, injured a security officer. Local firefighters and police responded, and the injured officer was treated and released from the hospital. On Tuesday, Captain Paul Starks, a local police spokesperson, indicated the chemicals from the explosion were consistent with those used for manufacturing methamphetamine. However, Starks did not indicate whether any equipment used in meth labs was found or whether NIST normally stores the chemicals used on site.
July 22, 2015 11:30 AM

### HP plans to power data centers with renewable energy

New York Times: On Tuesday, Hewlett-Packard became the latest technology giant to announce a plan to purchase power from renewable energy sources to supply major facilities. HP plans to buy 112 MW from a wind farm in Texas to supply its data centers in the state. Gabi Zedlmayer, HP's chief progress officer, said that the deal allows the company to reach its goal of reducing its greenhouse gas emissions by 20% from 2010 levels by 2020. SunEdison, the company that is providing the energy, will use the deal with HP to begin construction on another wind farm that will provide 300 MW when it is completed.

July 22, 2015 11:00 AM

### Dawn observes haze on Ceres

Nature: A mysterious haze has been detected hovering over the Occator crater on the dwarf planet Ceres. The observation was made by NASA’s Dawn spacecraft, which in March detected bright spots peppering Ceres’s surface. Occator has one large bright spot at its center and several smaller ones scattered nearby. Mission scientists have proposed that the bright spots are composed of some highly reflective material, such as ice or salt; the presence of the haze seems to point to sublimating ice. It is hoped that as Dawn’s IR spectrometer continues to map Ceres’s surface, it will eventually be able to determine definitively what it is.
July 21, 2015 3:30 PM

### Climate change may be accelerating more rapidly than predicted

Washington Post: A new study by climatologist James Hansen predicts that global warming will have more severe effects more quickly than previously thought. Along with 16 other researchers, Hansen says there is an “amplifying feedback” in relation to polar ice loss: As ice melts in Antarctica and Greenland, large volumes of cold, fresh water get dumped into the oceans, which slows down circulation patterns and causes a larger temperature contrast between the tropics and the poles. Increased ice melt from warm water trapped beneath the polar regions and more powerful storms created by the ever-increasing temperature differential will lead to greater sea-level rise and more catastrophic weather events than predicted by the Intergovernmental Panel on Climate Change in its 2013 climate assessment. That report had assumed a more linear process. The Washington Post article features feedback on the controversial study from a number of top climate researchers.
July 21, 2015 2:00 PM

### Study looks at problems experienced in robotic surgery

MIT Technology Review: According to a recent study, most of the robotic surgical procedures performed over the past 14 years have gone smoothly. However, a significant number have suffered some sort of adverse event, even if it did not result in the injury or death of the patient. Jai Raman at the Rush University Medical Center in Chicago and his colleagues analyzed data on robotic procedures recorded in the Manufacturer and User Facility Device Experience, or MAUDE, database kept by the US Food and Drug Administration. They found that of the more than 1 million robotic procedures carried out in the US between 2000 and 2013, 144 people died. However, procedure complications and prolonged operating times were caused by device malfunctions, which included equipment sparking and burning patients, pieces breaking off and falling into the patient’s body, equipment moving uncontrollably, and video feed being lost. Although no comparison was made with procedures that don’t use robotic techniques, the researchers say there is room for improvement in robotic equipment design and in the way accidents are reported and investigated.
July 21, 2015 11:45 AM

### SpaceX isolates possible cause of Falcon 9 rocket explosion

New York Times: On 28 June a Falcon 9 rocket taking supplies to the International Space Station (ISS) exploded about 139 seconds after liftoff. SpaceX CEO Elon Musk says the explosion was likely caused by the failure of a metal strut holding down a helium bottle in the upper-stage liquid oxygen tank. Although the metal struts had been certified to withstand up to 10 000 pounds of force, the one in question appears to have failed at just 2000 pounds. SpaceX has announced it will cease to use those particular struts and will perform additional hardware quality checks. A week after the SpaceX rocket failed, a Russian resupply ship succeeded in reaching the three astronauts on board the ISS. Before the 28 June launch failure, SpaceX’s Falcon 9 had successfully completed 18 missions.
July 21, 2015 10:50 AM

BBC: Not only are fossil-fuel emissions affecting Earth’s climate, but now it has been shown that they may diminish the effectiveness of an important technique for determining the age of objects—radiocarbon dating. Radiocarbon, or carbon-14, is created in the atmosphere when cosmic rays interact with atmospheric nitrogen. It is then absorbed by plants via photosynthesis and by animals that eat the plants. When the plant or animal dies, it stops absorbing 14C and the remaining 14C begins to decay. That decay rate can be used to calculate the organism’s age. Radiocarbon dating was developed in the late 1940s and has been used for various purposes, including the dating of ancient artifacts. With so much nonradioactive carbon being pumped into the atmosphere, however, the ratio of radioactive to nonradioactive carbon is becoming diluted, which makes objects appear to be older than they really are. That artificial aging of the atmosphere, which is worsening as fossil-fuel emissions increase, could have implications for a range of fields, including Earth science, archaeology, forgery detection, isotope forensics, and physiology, according to a study published in the Proceedings of the National Academy of Sciences.
July 20, 2015 2:42 PM

### Cyclical drying could have led to formation of earliest amino acids

Ars Technica: The process that led to the formation of the complex molecules and the development of life on Earth is still unknown. One of the big questions involves the formation of proteins, which are long chains of amino acids. Amino acids combine an acidic molecule with a nitrogen-based amine molecule. Two amino acids can bond to each other in a stable molecule, but the reaction is not energetically favorable because it produces an excess water molecule. However, increasing the concentration of amino acids in a solution increases the possibility of the reaction occurring spontaneously. A similar process happens with another molecule, lactic acid, but at a more energetically favorable rate. A new paper suggests that the drying of a solution that includes both amino acids and lactic acids could result in the two different molecules bonding together. However, over repeated cycles of drying and wetting, the amino acids replace lactic acids in the molecular chains. When the researchers tested the process, after just 20 cycles, 65% of the chains in the solution were made purely of amino acids.
July 20, 2015 2:10 PM

### Entrepreneur collaborates with scientists to search for alien life

BBC: A large-scale scientific search for signs of intelligent life in the universe is being launched by a private organization, Breakthrough Initiatives, founded by Russian entrepreneur Yuri Milner, who has a degree in theoretical physics. The group claims the \$100 million effort will be the largest such venture to date. It will make use of two of the world’s most powerful telescopes—the Green Bank Telescope in West Virginia and the Parkes Telescope in New South Wales, Australia—to scan 10 times more of the sky than previous projects and 5 times more of the radio spectrum. Among those involved with the initiative are renowned physicist Stephen Hawking and UK Astronomer Royal Martin Rees. "The search for extra-terrestrial life is the most exciting quest in 21st-century science,” said Rees.
See more
/content/aip/magazine/physicstoday/news/news-picks
201507
This is a required field


text:
# Training Neural Networks with Theano

Training neural networks involves quite a few tricky bits. We try to make everything clear and easy to understand, to get you training your neural networks as quickly as possible.

Theano allows us to write relatively concise code that follows the structure of the underlying maths. To run the code yourself, download the notebook at https://github.com/ASIDataScience/training-neural-networks-notebook

## Recognising hand-written digits

We will train a network to classify digits. More precisely, we want a network that when presented with an image of a hand-written digit will tell us what digit it is ('0', '1', ..., '9'). The data we will use for this task is known as the MNIST dataset. It has a long tradition in neural networks research, as the dataset is quite small but still very tricky to classify correctly.

### Import Theano

You can install theano with the Python package manager pip. At the command line type

pip install theano

or check out the theano documentation if you run into trouble

The theano.tensor module contains useful functionality for manipulating vectors and matrices, like we will be doing here, so let's import it along with the full package.

import theano
import theano.tensor as T
import numpy as np

With the data loaded into a variable mnist_dataset, we split it into three part: a training set, used to teach the network to recognize digits, a validation set that could be used to tune and compare models, and finally a test set, used to see how well it has learned.

We then prepare the datasets, splitting each set into the images and their labels, and store them in Theano shared variables, a bit of theano magic that's explained later. Understanding this massaging of the data isn't crucial.

train_set, valid_set, test_set = mnist_dataset
prepare_data = lambda x: (theano.shared(x[0].astype('float64')),
theano.shared(x[1].astype('int32')))
(training_x, training_y), (valid_x, valid_y), (test_x, test_y) = map(prepare_data, (train_set, valid_set, test_set))

Let's look at the first three images in the training set, then set about building a machine learning model that has the ability to recognise digits itself!

We've defined a function plot_mnist_digit exactly for printing out the training images - it's just for prettiness.

for i in range(3):
plot_mnist_digit(training_x.get_value()[i])             # training_x is a theano object containing *images*
print ('Image class: '+ str(training_y.get_value()[i])) # training_y is a theano object containing *labels*
Image class: 5
Image class: 0
Image class: 4

# The neural network model

The neural network model that we will build defines a probability distribution

P(Y = y | X = x; θ),

where Y represents the image class, which means it is a random variable that can take the values 0-9, X represents the image pixels and is a vector-valued random variable (we collapse the image matrix into a vector), and θ is the set of model parameters that we are going to learn.

In this tutorial we build two models: first implementing a logistic regression model, then extending it to a neural network model.

## Multi-class logistic regression

The equation for our first model is give by

P(Y = y | x; θ) ∝ [σ(xTW+bT)]y,

where [x]i is the ith entry of vector x, and the use of the proportionality symbol means that the probability is equal to the expression on the right hand side times a constant chosen such that yP(Y = y | x; θ) = 1

The parameter set θ for this model is θ = {W, b}, where W is a matrix and b is a vector. We also use the non-linearity σ given by

$$\sigma(t) = \frac{1}{1+e^{-t}}$$

Image credit: https://en.wikipedia.org/wiki/Logistic_function

When applied to a vector or matrix, the sigmoid function σ(t) is applied entrywise.

## Understanding the logistic regression model

One way to think about the logistic regression model is that it takes the input (x), puts it through a linear combination (xTW + bT) and then finally through a non-linearity: σ(xTW + bT).

The result of this operation is a vector representing the entire discrete distribution over all the possible classes - in our case the ten possible digits 0-9. To get the probability of a particular class y = 6 we extract the 6th entry of the probability vector: [σ(xTW+bT)]y

Graphically the model looks like this:

Each indiviual entry of the vectors x and P(Y) is shown as a circle -- known as the units (or artificial neurons) of the network. We have D input units (the dimensionality of the input vector, which is the flattened matrix of pixels) and C output units (the number of classes, which is the digits 0-9). The model parameters W and b are represented as arrows. We also show the application of the sigmoid functions, but we do not represent the normalization that makes the probabilities sum up to 1.

Another way to write the model above is using the SoftMax function. A good exercise is deriving the SoftMax function based on the fact that we can also write the same model using SoftMax and the equal sign instead of the proportionality sign:

P(Y = y | x; θ)=[SoftMax(Wx+b)]y,

## Neural network with a single hidden layer

Our second model is given by

$$P(Y = y \ |\ \boldsymbol{\textrm{x}} ; \theta) = \left[ \textrm{SoftMax} \left( \boldsymbol{\textrm{h}} ^\boldsymbol{\textrm{T}} \boldsymbol{\textrm{W}}_\boldsymbol{\textrm{hy}} + \boldsymbol{\textrm{b}}^\boldsymbol{\textrm{T}} _\boldsymbol{\textrm{hy}}\right)\right]_y, \\ \boldsymbol{\textrm{h}} = \tanh \left( \boldsymbol{\textrm{x}} ^\boldsymbol{\textrm{T}} \boldsymbol{\textrm{W}}_\boldsymbol{\textrm{xh}} + \boldsymbol{\textrm{b}}^\boldsymbol{\textrm{T}}_\boldsymbol{\textrm{xh}}\right)$$

This is very similar to the logistic regression model. Here we have introduced a new vector-valued random variable h. We call this a 'hidden' variable or 'latent' variable, as we do not have any data observations for it. This variable may not even correspond to any quantity in the real world, but we use it to increase the power of our statistical model.

We now also have more parameters: θ = {Wxh, Why, bxh, bhy}.

tanh is the hyperbolic tangent function given by
$$\tanh(t) = \frac{e^t-e^{-t}}{e^t+e^{-t}}$$

Image credit: http://mathworld.wolfram.com/HyperbolicTangent.html

Like with the sigmoid, when applied to a vector or matrix, tanh function is applied entrywise.

Graphically, this model looks like this:

Now our depiction shows a hidden layer with M units (this number can be different from the number of input neurons and number of output neurons), and we have two different nonlinearities in the graph: tanh and sigmoids (but again we are not graphically representing the SoftMax normalization).

## Teaching the network

The way we're going to make our network learn is by trying to find some values for our parameters θ so that the network is as likely as possible to guess the correct class. That is, given a data set of training images x1, x2, …, xN and correct labels y1, y2, …, yN, we want to find the parameters that maximize probability of the correct labels given the images.

This method of choosing parameters is called maximum likelihood (ML), and we can express it mathematically as finding the parameters θ which maximize the likelihood function:

θ* = argmaxθ P(Y1 = y1, Y2 = y2, …, YN = yN | X1 = x1, X2 = x2, …, XN = xN; θ)

And since our data points are independent, we can write this joint probability as a product of probabilities.

## Writing the likelihood function for an entire dataset

In our likelihood function above, each random variable pair (X1, Y1),(X2, Y2),…,(XN, YN) refers to a single data point. But since virtually all of our computations need to deal with multiple data points, we will find it both useful and computationally efficient to express both the mathematics and our Python code in terms of datasets. Thus, we will express the scalar random variables Y1, Y2, …, YN as the vector-valued random variable Y, and the vector-valued random variables X1, X2, …, XN as the matrix-valued random variable X, where the matrix X has as many rows as there are data points.

Using this notation, we rewrite the maximum likelihood equation above:

θ* = argmaxθ P(Y = y | X = X; θ)

Similarly, we can specify the logistic regression model it terms of multiple datapoints:

P(Y | X = X; θ)=SoftMax(XW+1bT)

Here the result is a matrix of probabilities with as many rows as there are data points, and as many columns as there are classes. We also consider the SoftMax to normalize the result of the linear combination (XW+1bT) in such a way that each row of the result is a proper probability distribution summing to 1. Note that we have had to multiply the bias vector b by the vertical vector of all ones 1 in order to add the bias term for every single data point.

The neural network model equations follow a similar pattern, which it would be a good exercise to write out for yourself.

## Computing the log-likelihood of a dataset

In most machine learning applications, it is better to maximize the log-likelihood rather than the likelihood. This is done because the log-likelihood tends to be simpler to compute and more numerically stable than the likelihood. In terms of the math, this doesn't make things much more complicated, as all we need to add is a log in front of the likelihood:

θ* = argmaxθ logP(Y = y | X = X; θ)

Since the logarithm of a product is the sum of the logarithms, we can write:

$$\log P( Y = \boldsymbol{\textrm{y}}\ | \ X = \boldsymbol{\textrm{X}} ; \theta) = \sum_{i=1}^N \log P( Y_i = \boldsymbol{\textrm{y}}_i\ | \ X_i = \boldsymbol{\textrm{X}}_i ; \theta)$$

that is, the log joint probability is the sum of the log marginal probabilities.

Now let's plug in the probability of a dataset from above to obtain:

logP(Y = y | X = X; θ)=∑[log(SoftMax(XW+1bT))]⋅, y

where we use the notation [M]⋅, y to mean that from the matrix M we construct a new vector (a1, a2, …, an) such that ai = Mi, yii. We use a slight abuse of notation and we use the sum symbol to indicate the summation of the entries in the vector; we need a summation because the log joint probability is the sum of the log marginal probabilities.

## Classifying a dataset

Once we have a set of parameters θ that we are happy with, we can use the model to classify new data.

We have built a model that gives us a distribution over classes given the data. How do we assign a class to a new data point? The simplest way is to choose the class with the highest probability under the model to be the class we assign. We can write this mathematically, again using vector notation:

$$\hat{\boldsymbol{\textrm{y}}} = {\arg\max} \ P(Y \ |\ \boldsymbol{\textrm{X}} ; \theta)$$

Computing the maximum likelihood parameters is computationally unfeasible, so we're going to use a method called gradient ascent to find a set of parameters that are really good but perhaps not the absolute best.

The idea of gradient ascent is very simple. Given a function f : ℝn → ℝ, we want to iteratively find points f(xn) such that the value of the function gets progressively higher, that is: f(xn + 1)>f(xn)∀n. One way to do this is taking the direction of steepest ascent, which is just the gradient of the function $\frac{\partial f}{\partial \boldsymbol{\textrm{x}}}$, and taking a step in that direction times a constant λ known as the learning rate that describes how big the step should be. We express this mathematically as:

$$\boldsymbol{\textrm{x}}_{n+1} \leftarrow \boldsymbol{\textrm{x}}_n + \lambda \frac{\partial f}{\partial \boldsymbol{\textrm{x}}}$$

The last detail is choosing the starting point, x0, which we can arbitrarily choose by setting to zero or to some random value. Graphically, the algorithm looks like this, with each color representing the path from a different starting point:

Note that this method tends to find the top of the nearest hill ('local' maximum), and not the overall best point ('global' maximum).

It is also not guaranteed to increase the value of the function at each step; if the learning rate is too large, the algorithm could potentially jump across the top of the hill to a lower point on the other side of the hill. Much research goes into optimization methods, and many neural networks models are trained with methods that are more complicated than gradient ascent as it's presented here, but this same idea is at the base of all of those methods.

Finally, most people talk about and use gradient descent on the negative log-likelihood rather than gradient ascent on the log-likelihood; this is because gradient descent is the standard algorithm in the field of optimization. Gradient ascent in used in this tutorial to keep things a bit simpler.

## Optimizing our model with gradient ascent

This is extremely easy: we just apply the equation above to our parameters taking the gradient of the log-likelihood:

\begin{align} \boldsymbol{\textrm{W}} &\leftarrow \boldsymbol{\textrm{W}} + \lambda \frac{\partial \log P( Y = \boldsymbol{\textrm{y}}\ | \ X = \boldsymbol{\textrm{X}} ; \theta)}{\partial \boldsymbol{\textrm{W}}} \\ \boldsymbol{\textrm{b}} &\leftarrow \boldsymbol{\textrm{b}} + \lambda \frac{\partial \log P( Y = \boldsymbol{\textrm{y}}\ | \ X = \boldsymbol{\textrm{X}} ; \theta)}{\partial \boldsymbol{\textrm{b}}} \\ \end{align}

# Coding with Theano

Coding with Theano is extremely simple: once we have the equations behind the model, we pretty much type them directly. Since we will be training on multiple data points, we are going to encode the model as we wrote it for datasets, using vectors and matrices.

P(Y | X; θ)=SoftMax(XW+1bT)

The first thing we do is declare all of the variables (X, y, W, b) that we will be using and their types like you would in Java or C:

n_classes = 10  # each digit is one of 0-9
dims = 28 * 28  # our input data is flattened 28x28 matrices of image pixels

X = T.dmatrix() # Theano double matrix
y = T.ivector() # Theano integer vector
W = theano.shared(np.zeros([dims,n_classes])) # Theano shared double matrix
b = theano.shared(np.zeros(n_classes))        # Theano shared double vector

As you can see above, we defined W and b to be shared variables. This means that the values of these variables are persistent -- their values live on after we have run Theano operations. This is opposed to regular Theano variables, which only take values when Theano runs, and otherwise only exist in the abstract.

The reason for making W and b shared variables is that we want to run multiple iterations of gradient descent, and to do that, we need their values to persist. Furthermore, we want to find good parameters through training, but we will then want to use the same parameters for prediction, so we need them to be persistent

Let's now write our statistical model in Theano. We basically copy the following equation into code:

P(Y | X; θ)=SoftMax(XW+1bT)

P = T.nnet.softmax(T.dot(X,W) + b) # the matrix of probabilities of all classes for all data points

Theano provides us with a T.nnet.softmax function to compute SoftMax, correctly normalizing the probabilities so that each row of the matrix P is a proper probability distribution that sums to 1.

Note that we didn't need to multiply b by the 1 vector, Theano will do the correct addition for us automatically, just like numpy would do it. Here is a simple numpy example illustrating this:

amatrix = np.zeros((3,2)) # 3x2 matrix of all zeros
avector = np.array((1,2)) # the vector [1,2]
amatrix + avector
array([[ 1.,  2.],
[ 1.,  2.],
[ 1.,  2.]])

Our next equation is the log-likelihood (LL) of a dataset:

logP(Y = y | X = X; θ)=∑[log(SoftMax(XW+1bT))]⋅, y

LL = T.mean(T.log(P)[T.arange(P.shape[0]), y]) # the log-likelihood (LL)

OK... there's a lot going on here.

We have used two important tricks: - using mean instead of sum, - using the strange indexing [T.arange(P.shape[0]), y]

Let's go over them one by one.

#### Use the mean instead of the sum

Imagine that we were to construct a new dataset that contains each data point in the original dataset twice. Then the log-likelihood of the new dataset will be double that of the original. More importantly, the gradient of the log-likelihood of the new dataset will also be double the gradient of the log-likelihood of the original dataset. But we would like the size of the gradient to not depend on the amount of duplication in our dataset, and the easiest way to accomplish that is to divide the gradient by the size of the dataset.

Since taking the mean is equivalent to taking the sum and then dividing by the number of data points, what we are computing here is a type of "normalized" log-likelihood that will cause our gradient descent algorithm to be robust to change in dataset size. The quantity we are computing in the code can be more precisely described as the average log-likelihood for a single datapoint.

#### Use indexing to create a vector from a matrix

The second thing we do is use [T.arange(P.shape[0]), y] to apply the mathematical operation denoted in the equation by [⋅]⋅, y, that is, constructing a new vector (a1, a2, …, an) from a matrix M such that ai = Mi, yii.

As cryptic as it may be, this is a peculiar, but standard numpy way to index. For example, given a matrix

**
M = np.array([[1,2,3,4],
[5,6,7,8],
[9,10,11,12]]

**

If we wanted to extract one element from each row, say the 1st element from the first row, and the last from the others

**
M[0,0], M[1,3], M[2,3]

**

We could write that as a single index expression by combining the indexes

**
M[(0,1,2), (0,3,3)]

**

But now the first index is just (0…#rows − 1), or, in code, np.arange(M.shape[0]).

So we can write:

M = np.array([[1,  2,  3,  4],
[5,  6,  7,  8],
[9, 10, 11, 12]])

M[np.arange(M.shape[0]), (0,3,3)]
array([ 1,  8, 12])

We're done with the model. There's one more thing we need to do, and that is specify the gradient updates

learning_rate = 0.5 # we tuned this parameter by hand
[W, W + learning_rate * T.grad(LL, W)],
[b, b + learning_rate * T.grad(LL, b)]
]

## Theano functions

OK, we're done coding the model. What do we do next?

When working with Theano, the next step is to create a Theano function. A Theano function is the basic unit of Theano code that we call to do something. In our case, this something will be performing a single gradient ascent iteration.

We create a Theano function by calling the function theano.function (yes, we create a function by calling a function). theano.function has four important parameters that we provide in order to get a Theano function:

• inputs -- the list of input variables of the Theano function, similar to the inputs of a Python function
• outputs -- the list of output variables of the Theano function, similar to the inputs of a Python function
• updates -- a list of updates for shared variables, in the format we used above when we defined the variable updates
• givens -- a dictionary that allows substituting some variables from the model with other variables

In our case, we want the input to be the training dataset, the updates to be the gradient ascent updates, and while we don't really need an output, it will be helpful to get the log-likelihood as an output to see that we are doing the right things.

However, we will use the givens instead of the input to provide the data to the function.

Doing it this way is more efficient, as we've already loaded up the training dataset into memory as a shared Theano function, when we first loaded the data.

Our Theano function will look like this:

training_function = theano.function(
inputs = [],             # use givens instead of the inputs as it's more efficient
outputs = LL,            # output log-likelihood just to check that it is improving
givens = {X: training_x, # we indicate that the model variables X and y defined in the abstract
y: training_y} # should take the values in the shared variables training_x and training_y
)

OK, let's run ten iterations of our code and see what the log-likelihood does

for i in range(10):
current_LL = training_function()
print("Log-likelihood = " + str(current_LL) + "\t\t" +
"Average probability of the correct class = " + str(np.exp(current_LL))
)
Log-likelihood = -2.302585093       Average probability of the correct class = 0.0999999999998
Log-likelihood = -1.83912090251     Average probability of the correct class = 0.158957103494
Log-likelihood = -1.52505436796     Average probability of the correct class = 0.217609225573
Log-likelihood = -1.31424405922     Average probability of the correct class = 0.268677350658
Log-likelihood = -1.16824101592     Average probability of the correct class = 0.310913352197
Log-likelihood = -1.06176739247     Average probability of the correct class = 0.34584402773
Log-likelihood = -0.98204553767     Average probability of the correct class = 0.374544170518
Log-likelihood = -0.919423093595        Average probability of the correct class = 0.398749015602
Log-likelihood = -0.869394886438        Average probability of the correct class = 0.419205139229
Log-likelihood = -0.827688588441        Average probability of the correct class = 0.437058341404

## Using the model for classification

Great, it appears we're improving the log-likelihood of the model on the training data. Now to use it to classify. Recall our equation that expresses how to use the model to get the class:

$$\hat{\boldsymbol{\textrm{y}}} = {\arg\max} \ P(Y \ |\ \boldsymbol{\textrm{X}} ; \theta)$$

Let's put that in code:

y_hat = T.argmax(P, axis=1)

Note that we had to specify axis=1, that is, we want to get the argmax for each row, as each row represents the distribution for one datapoint.

### Create a Theano function that classifies a dataset

Similarly to the training function, the classification function will use givens to pass in the test dataset, and output y_hat which we just defined

classification_function = theano.function(
inputs = [],
outputs = y_hat,
givens = {X:test_x} # don't need the true labels test_y here
)

Now let's run the classification once, and print the first ten images, the true labels, and the labels assigned by the model

test_y_hat = classification_function()
print ("Classification error: "+ str(100 * (1 - np.mean(test_y_hat  test_y.get_value()))) + "%")
for i in range(10):
plot_mnist_digit(
test_x.get_value()[i]                                # test_x is a theano object of images
)
print ('Image class: \t\t' + str(test_y.get_value()[i])) # test_y is a theano object of *true labels*
print ('Model-assigned class: \t' + str(test_y_hat[i]))  # test_y_hat is a theano object of *predicted labels*
Classification error: 15.12%
Image class:        7
Model-assigned class:   7

Image class:        2
Model-assigned class:   2
Image class:        1
Model-assigned class:   1
Image class:        0
Model-assigned class:   0
Image class:        4
Model-assigned class:   4
Image class:        1
Model-assigned class:   1
Image class:        4
Model-assigned class:   4
Image class:        9
Model-assigned class:   9
Image class:        5
Model-assigned class:   2
Image class:        9
Model-assigned class:   9

## Training a neural network

So far we have trained a logistic regression model. The neural network model is so similar that we can imlepement it with just a few changes to the code.

We need
- to declare the hidden layer variable - to decide on the size of the hidden layer (we'll keep it small so it will run on your personal computer) - new parameters

H = T.dmatrix()  # Theano double matrix
hidden\_layer\_size = 20
W_xh = theano.shared(0.01 * np.random.randn(dims, hidden_layer_size))
W_hy = theano.shared(np.zeros([hidden_layer_size, n_classes]))
b_xh = theano.shared(np.zeros(hidden_layer_size))
b_hy = theano.shared(np.zeros(n_classes))

Remember our model? Let's write it using matrices and vector for an entire dataset:

$$\boldsymbol{\textrm{H}} = \tanh \left( \boldsymbol{\textrm{X}} \boldsymbol{\textrm{W}}_\boldsymbol{\textrm{xh}} + \boldsymbol{\textrm{1b}}^\boldsymbol{\textrm{T}}_\boldsymbol{\textrm{xh}}\right), \\ P(Y \ |\ \boldsymbol{\textrm{x}} ; \theta) = \textrm{SoftMax} \left( \boldsymbol{\textrm{H}} \boldsymbol{\textrm{W}}_\boldsymbol{\textrm{hy}} + \boldsymbol{\textrm{1b}}^\boldsymbol{\textrm{T}} _\boldsymbol{\textrm{hy}}\right)$$

Let's code it up:

H =         T.tanh( T.dot(X, W_xh) + b_xh )
P = T.nnet.softmax( T.dot(H, W_hy) + b_hy )

LL = T.mean(T.log(P)[T.arange(P.shape[0]), y])  # the log-likelihood (LL)

[W_xh, W_xh + learning_rate * T.grad(LL, W_xh)],
[W_hy, W_hy + learning_rate * T.grad(LL, W_hy)],
[b_xh, b_xh + learning_rate * T.grad(LL, b_xh)],
[b_hy, b_hy + learning_rate * T.grad(LL, b_hy)],
]

### Redefining Log-loss

The one extremely important thing we did here was redefined LL.

This is a crucial point about how Theano works:

Whenever we define a theano variable, like we did with P, we create a new object. When we define a new theano variable in terms of another Theano variable, like we did with LL, using **
LL = T.mean(T.log(P)[T.arange(P.shape[0]), y])

we implicitly create a new object for LL** that has a reference to our variable P we just defined.

Now the crucial part: say we redefine P. Then our variable LL still has a reference to the old variable P, and we need to update the reference to LL by re-running the definition for LL for everything to work correctly.

Bugs in Theano are very commonly produced by exactly this. It is a good reason to always use Theano in scripts rather than in a notebook like we are here.

Phew, we are now ready to train!

training_function = theano.function(
inputs = [],
outputs = LL,
givens = {X: training_x[:5000],  # use only 10% of the data so model not too complicated
y: training_y[:5000]}  # to train on a personal computer
)
for i in range(60):                  # train more than for logistic regression as this model is more complex
current_LL = training_function()
print(
"Log-likelihood = " + str(current_LL) + "\t\t" +
"Average probability of the correct class = " + str(np.exp(current_LL))
)
Log-likelihood = -2.30258509299     Average probability of the correct class = 0.1
Log-likelihood = -2.30123532532     Average probability of the correct class = 0.100135067902
Log-likelihood = -2.2998167966      Average probability of the correct class = 0.100277213167
Log-likelihood = -2.29812291581     Average probability of the correct class = 0.100447214752
Log-likelihood = -2.29590367206     Average probability of the correct class = 0.100670379142
Log-likelihood = -2.29281799209     Average probability of the correct class = 0.100981495471
Log-likelihood = -2.28837617693     Average probability of the correct class = 0.101431034253
Log-likelihood = -2.28186920166     Average probability of the correct class = 0.102093195481
Log-likelihood = -2.27228942655     Average probability of the correct class = 0.103075924982
Log-likelihood = -2.25826075997     Average probability of the correct class = 0.104532133214
Log-likelihood = -2.23801271431     Average probability of the correct class = 0.1066702782
Log-likelihood = -2.20944653423     Average probability of the correct class = 0.109761380875
Log-likelihood = -2.17036118024     Average probability of the correct class = 0.114136385658
Log-likelihood = -2.11893562814     Average probability of the correct class = 0.120159454814
Log-likelihood = -2.05448392668     Average probability of the correct class = 0.128158957933
Log-likelihood = -1.9781664624      Average probability of the correct class = 0.138322624675
Log-likelihood = -1.89308937088     Average probability of the correct class = 0.150605812178
Log-likelihood = -1.80356444392     Average probability of the correct class = 0.16471073844
Log-likelihood = -1.71390604967     Average probability of the correct class = 0.180160699808
Log-likelihood = -1.62743071587     Average probability of the correct class = 0.196433620113
Log-likelihood = -1.54612544874     Average probability of the correct class = 0.213071934689
Log-likelihood = -1.47089832383     Average probability of the correct class = 0.22971903039
Log-likelihood = -1.40197788868     Average probability of the correct class = 0.246109704628
Log-likelihood = -1.33918043475     Average probability of the correct class = 0.262060356155
Log-likelihood = -1.28205618089     Average probability of the correct class = 0.27746619282
Log-likelihood = -1.23000534858     Average probability of the correct class = 0.292291014336
Log-likelihood = -1.1823856608      Average probability of the correct class = 0.306546549485
Log-likelihood = -1.13858981262     Average probability of the correct class = 0.320270344716
Log-likelihood = -1.09808364787     Average probability of the correct class = 0.333509593518
Log-likelihood = -1.06041545036     Average probability of the correct class = 0.346311905034
Log-likelihood = -1.02521125894     Average probability of the correct class = 0.358720674451
Log-likelihood = -0.992165825303        Average probability of the correct class = 0.37077279169
Log-likelihood = -0.96103293577     Average probability of the correct class = 0.382497586412
Log-likelihood = -0.931615847613        Average probability of the correct class = 0.393916686506
Log-likelihood = -0.903757946709        Average probability of the correct class = 0.405044659855
Log-likelihood = -0.877333927776        Average probability of the correct class = 0.415890228318
Log-likelihood = -0.852241954255        Average probability of the correct class = 0.426457760591
Log-likelihood = -0.828397170186        Average probability of the correct class = 0.436748759536
Log-likelihood = -0.805726712165        Average probability of the correct class = 0.446763140352
Log-likelihood = -0.784166136418        Average probability of the correct class = 0.456500202015
Log-likelihood = -0.76365701148     Average probability of the correct class = 0.465959288932
Log-likelihood = -0.744145358272        Average probability of the correct class = 0.475140201108
Log-likelihood = -0.725580637598        Average probability of the correct class = 0.484043433528
Log-likelihood = -0.707915058331        Average probability of the correct class = 0.492670316262
Log-likelihood = -0.691103068952        Average probability of the correct class = 0.501023101114
Log-likelihood = -0.675100970071        Average probability of the correct class = 0.509105013643
Log-likelihood = -0.65986663098     Average probability of the correct class = 0.516920271045
Log-likelihood = -0.645359309369        Average probability of the correct class = 0.524474059803
Log-likelihood = -0.631539569966        Average probability of the correct class = 0.531772469537
Log-likelihood = -0.61836928705     Average probability of the correct class = 0.538822386202
Log-likelihood = -0.605811706762        Average probability of the correct class = 0.545631354182
Log-likelihood = -0.593831541869        Average probability of the correct class = 0.552207420302
Log-likelihood = -0.582395074029        Average probability of the correct class = 0.558558973142
Log-likelihood = -0.571470244429        Average probability of the correct class = 0.564694589
Log-likelihood = -0.561026720601        Average probability of the correct class = 0.570622892704
Log-likelihood = -0.551035933447        Average probability of the correct class = 0.576352438247
Log-likelihood = -0.541471083378        Average probability of the correct class = 0.581891611356
Log-likelihood = -0.532307117743        Average probability of the correct class = 0.587248554016
Log-likelihood = -0.523520683596        Average probability of the correct class = 0.592431109513
Log-likelihood = -0.515090060711        Average probability of the correct class = 0.597446785713

## Test the model

y_hat = T.argmax(P, axis=1)
test\_y\_hat = classification_function()

print ("Classification error: " + str(100 * (1 - np.mean(test_y_hat  test_y.get_value()))) + "%")

for i in range(10):
plot_mnist_digit(
test_x.get_value()[i]                                # test_y is a theano object of *images*
)
print ('Image class: \t\t' + str(test_y.get_value()[i])) # test_y is a theano object of *true labels*
print ('Model-assigned class: \t' + str(test_y_hat[i]))  # test_y_hat is a theano object of *predicted labels*
Classification error: 15.12%
Image class:        7
Model-assigned class:   7
Image class:        2
Model-assigned class:   2
Image class:        1
Model-assigned class:   1
Image class:        0
Model-assigned class:   0
Image class:        4
Model-assigned class:   4
Image class:        1
Model-assigned class:   1
Image class:        4
Model-assigned class:   4
Image class:        9
Model-assigned class:   9
Image class:        5
Model-assigned class:   2
Image class:        9
Model-assigned class:   9

## We've trained a neural network with Theano!

Along the way we used a good chunk of mathematics, a hand-full of tricks but required very few lines of code.

Hopefully this leads nicely to using Theano more, in particular we would be excited to see you set theano to work with: smarter training algorithms, parallelized training, and training more popular network architectures.

Good luck!

By Gabi Teodoru.


text:
# Non-Stationary: Larger-than-unit root [duplicate]

I keep reading everywhere that a time series is non-stationary (e.g. http://en.wikipedia.org/wiki/Unit_root or http://en.wikipedia.org/wiki/Stationary_process) if there's a unit-root. But isn't a larger-than-one root also means the series is non-stationary? Why don't people say that roots larger than unit imply non-stationarity instead?

Intuitive explanation of unit root is a much more general question and indeed the top answer has one sentence in the long essay addresses this question. For those who's only interested in the answer to my question and did not want to read the whole thing (which is epic by the way): "Economists are perhaps the greatest analysts of time series and employers of the AR process technology. Their series of data typically do not accelerate out of sight. They are concerned, therefore, only whether there is a characteristic direction whose value may be as large as 1 in size: a "unit root."

• You can have nonstationarity in vastly more circumstances than unit roots – Glen_b Mar 25 '15 at 11:20

I think this is actually a quite good question, which is often neglected (as you have noticed) and which I myself haven't thought about much before. The main point, I would say, is that processes with larger-than-one roots (called explosive roots) are not as interesting. If you have something which is just slightly above one, the process will fairly quickly just look like a nice curve. An explosive process will therefore reveal itself, but the (visual) difference between a unit root process and a near-unit root process is much more subtle.

Consider the AR(1) process $$y_t=ay_{t-1}+\epsilon_t.$$ I have simulated this with $a=1$ (this is the $y_t$ process in the figures), which is a random walk with a unit root. Also shown is $x_t$ which is the same as above but with a slight perturbation, so $a=1.05$ now. Thus, it has an explosive (not just a unit) root. As you can see, the behavior they exhibit is quite different (granted this is just one simulation, of course). You see the trending-like behavior already with $T=40$, and with $T=1000$ it just looks odd. Therefore, as I see it, you disregard the possibility of an explosive root many times because it is "unrealistic". A process such as what you have in the top right panel might instead, in practice, be modeled using deterministic trends with a possible non-stationary process moving around this trend.

So, non-stationarity is definitely implied by explosive roots. But in practice these are much less often found, so we spend quite some time learning about the more realistic situation of non-stationarity, which is a unit root. For the same reason, you often don't learn a whole lot about a negative unit root (i.e. $a=-1$ in the model above).

eps  <- rnorm(1000)
eps2 <- rnorm(1000)
y <- eps
x <- eps2
for (t in 2:1000) {
y[t] <- y[t-1] + eps[t]
x[t] <- 1.05*x[t-1] + eps2[t]
}

par(mfrow=c(2,2))
plot(y[1:40], type = "l", ylab = "y, t=1, ..., 40", main = "a = 1")
plot(x[1:40], type = "l", ylab = "x, t=1, ..., 40", main = "a = 1.05")
plot(y, type = "l", main = "a = 1")
plot(x, type = "l", main = "a = 1.05")

There are several kinds of non-stationarities:

1) Series expected value is a function of time
2) Series variance depends on time and not just about lag
3) Etc

Series with linear trend is non-stationary but stationary around trend..

EDIT:

Your example of stochastic difference equation with explosive root is of course non-stationary if you take definition from non time independent variance or expected value.

But linear trend is more interesting as a mathematical model than explosive stochastic difference equation.

• How exactly does this answer any of the two questions "But isn't a larger-than-one root also means the series is non-stationary?" or "Why don't people say that roots larger than unit imply non-stationarity instead?"? – hejseb Mar 25 '15 at 6:41
• @hejseb Added more on edit... – Analyst Mar 25 '15 at 7:11