text:
Forrest Moses '56 and "Sylvan Waters" in Santa Fe
By Julie Campbell
If your July travels find you in Santa Fe, N.M., you'll want to take in "Sylvan Waters," a new exhibition of paintings and monotypes by Forrest Moses, a member of Washington and Lee's Class of 1956.
As Forrest told the Albuquerque Journal in an article about the show, "I have loved water for years. I have been painting streams and ponds and lakes for a long time. I like the details of water places. Those details can be abstracted more easily than anything else, especially reflections. I prefer to work in reflections because that takes away the realism of the actual landscape."
On his detailed and handsome website, forrestmoses.com, you can view many examples of his work, which encompasses oil paintings, monotypes and photography. He's also combined his images with his words in a limited edition book, "Forrest Moses," printed in Italy and published in 2001.
Forrest was born in Virginia and received a B.A. from W&L. Following graduation, he served in the Navy, traveled in Europe and studied at the Pratt Institute, in New York City. He then lived and worked in Houston and on the Monterey Peninsula of California before settling in Santa Fe in 1969.
He's had more than 40 one-man shows around the country. This one, at Santa Fe's LewAllen Galleries, runs until July 29.


text:
Q: Manually add data to MySQL query array I've just started using Codeigniter, and am loving MVC. I have a simple mysql query in a model that returns an array of rows to my controller.
$query = $this->db->get('shows');
return $query->result();

Date information is stored in the database as mysql date (yyyy-mm-dd) and I explode to get month mm and day dd.
What I'm wondering is, is there any way to manually add the variables for month and day of each row to the query result using a simple foreach? I know I could manually add each database field's value to an array, and include the month and day variables, but I'm hoping there's a simpler way of inserting them into the already existing array created by the query.

A: Ok im not familiar with codeigniter, but if you want to get month and day in raw mysql you can simply do something like below
$sql = 'SELECT MONTH(date) AS month, DAY(date) AS day FROM shows';
$result = $this->db->query($sql);

Im assuming your date column is in date dataype

A: @slier had the code, but below is transformed to codeigniter
$this->db->select("*, MONTH(date_column) as date_month, DAY(date_column) as date_day");
$query = $this->db->get('shows');
return $query->result();

redpajama_set_name: RedPajamaStackExchange


text:
When Even the Sky Hurts
Arts and Culture, Poetry
Review by Mary Meriam
"The writing of a poem is also the act of taking a stand against the sadness."
[dropcap]M[/dropcap]y homeless sister once dreamily said that she would like a wardrobe, a closet full of nice clothes. Clothes are a home for our bodies. Clothes broadcast to others where we stand in the world. If our clothes are clean, comfortable, fashionable, and a good fit, we appear to be at home in the world.
Homelessness assaults a person from so many angles — hunger, thirst, cold, heat, danger, despair — that we sometimes can't help dreaming that if we only had the right clothes, we could have a home in the world. But while trying to hold together body and soul, our clothes get dirty, worn out, and torn, making us feel and look even more like misfits.
Jonathan Edwards (1703-1758) was an American preacher and philosopher who traveled alone on horseback from parish to parish. He had ideas and insights while he was traveling, and devised a method to remember his ideas by pinning small pieces of paper to his clothing. A contemporary poet, Susan Howe, was inspired by this story to write in her book, Souls of the Labadie Tract, "Words give clothing to hide our nakedness. I love to imagine this gaunt and solitary traveler covered in scraps riding through the woods and fields of Massachusetts and Connecticut."
Julia Vinograd's new book of poems, When Even the Sky Hurts, begins with "Shirt," which begins with these lines: "I'm wearing a dark cotton shirt / with yellow dancing skeletons, kicking high." Like the layers of clothes we put on to keep warm, this image has many layers of meaning. The shirt is worn by the speaker of the poem, whose words are both literal clothing for the speaker and figurative flesh for the skeletons. Skeletons emerge from a dark background to dance on a shirt; not just any skeletons, but cheerful "yellow" skeletons; and not just dancing, but "kicking high," full of energy.
Poet Lucia Perillo wrote that "the writing of a poem is also the act of taking a stand against the sadness." Vinograd's poems do more than take a stand; they dance away sadness, give us joyful clothing, transform death into life.
In the first few lines of "Letters to Soldiers," Vinograd compares "a new flowered dress" with "a mask of rain and blood" to show a heartbreaking contrast between two worlds that can no longer communicate:
Their girlfriends' letters, lightly perfumed
with gently looping penmanship
talk about a bake sale to support the troops
and describe a new flowered dress she bought
to match her strawberry shortcake
and he won't be such a great hero
that he forgets all about her, will he?
Soldiers read these letters, peering thru a mask
of rain and blood. Flies buzz over knife slashes.
Gentle, perfumed, flowered, strawberry words in the girlfriends' letters cannot give the soldiers protection from war's wounds and horrors — the soliders are wordless and exposed.
Clothing is transformed into food in these vivid, evocative lines from "Big Momma Blues":
"Baby, baby," sings Big Mama Blues, peeling
the city's orange vests off tired construction workers
like oranges fresh off a tree, peeling off soldiers' uniforms
like onion skins over a wound.
Big Mama's song brings comfort and relief, by stripping away the grown-up costumes of work and war, and making the workers and soldiers vulnerable, nourished babies, swaddled in song.
In the first two lines of a great poem, Vinograd answers the question "Why I Write Poetry": "Because I can't trust God / to look after the world and my friends." These words express Vinograd's wellspring of inspiration and compassion.
The middle of the poem continues with candor: "God, it is a very beautiful world, / but no thank you, it is not enough." The last words of the poem turn again toward the reader, with saving grace: "I'm throwing you a rope of words. / Hold on." Here words move beyond poems, letters, and songs, to the poet's naked voice providing a lifeline.
Even a celestial body has lost its clothes: "the moon's cast-off trousers." ("Mother"). Nakedness is cold, dark, impersonal, wordless: "in my own naked body. Night cold. Nameless." ("The Hunt").
Berkeley poet Julia Vinograd.
The following lines from "Street Crazy" describe a poignant relationship between the human and avian worlds, as if, when there's no other way to clothe our nakedness, we could find a home with birds, or wear the clothing of birds:
Crazy guy pulls swans out of the holes in his pocket.
They circle his head and fly away.
He picks up one long feather for his hat band
if he ever has a hat, but he drops the feather
a block later staring at a sale sign:
"everything must go." He feels himself going.
So sad, "if he ever has a hat," and he loses the feather he found, and everything is going. This guy isn't crazy to want nice clothes.
"Scorpion Dance" is reminiscent of the skeletons in "Shirt." Here a spider with a poisonous sting is alive on a t-shirt. Trouble is held at bay by the strong muscles of heartbeat, song, and loveliness.
The scorpion climbs over muscles of music,
bouncing against the banjo player's t-shirt
listening for his heartbeat in a storm of song,
in his lean body lovely between the claws of time.
I'd like to recommend Vinograd's poem, "First Love," in which there's nothing between the speaker and the wind that loves her. In "The Bogeyman's Love Song," clothes are no defense against the trouble of a seductive snake: "my long forked tongue will lick under your clothes." The woman, Jerusalem, tries to heal devastation and destruction in "Bridal Gowns and Jerusalem."
The most powerful and beautiful clothes are worn by this dancer, who makes the whole world new again. Here is the complete poem:
Bellydancer in a Coffeehouse
Not what I expected. Only her feet were bare
and her stilled face was naked.
At her slim waist deep blue strings spun
over a crimson whirling skirt with yellow satin harem pants
underneath. Brown lace sleeves shot with gold thread,
long and wide as trailing wings.
She started the dance kneeling under a scarf
and stood up so slowly the music
had to get under her knees and help.
Her body carved the air gradually out of her way.
Long dark hair piled in braids on her head
carried like a jug of water never spilling.
Her arms uncoiled but not to the crowd.
She never smiled.
The coffeehouse dissolved around her
it was too young and easy for her ruthless eyes.
She was an answer to genesis,
flesh become word.
Everything in Vinograd's world has feelings, even clothes, even skeletons. Earth is clothed by the sky, and even the sky hurts. I'm grateful for Vinograd's poems. Her words are like beautiful, comfortable clothing that she has given us to wear, so that we can feel at home in the world. You can find Julia Vinograd on Telegraph or at Fourth Street, hawking When Even the Sky Hurts.
Purchase When Even the Sky Hurts from Zietgeist Press in Berkeley from their website at http://ww.zeitgeist-press.com/vinograd.cfm

redpajama_set_name: RedPajamaCommonCrawl


David Cresswell does not yet have a rating. At this time, there are no reviews or comments for this chiropractor.
If you have personal experience with David Cresswell, we encourage you to share that experience with our ServiceRating.ca community. Your opinion is very important and David Cresswell will most certainly appreciate the feedback.


text:
Students seek more funding for illegal immigrant graduate students
Caitlin Grimes | North Carolina Campus Correspondent
Thursday, June 11, 2015 1:49 PM
Students at Claremont University, a graduate university in California, are pursuing more resources for graduate students residing in the United States illegally.
A recent story published in the LA Times highlighted the difficulties of student Iliana Perez in funding her graduate studies. A non-resident of the United States, Perez struggled to find funding for continuing her education until she received a scholarship from Claremont University.
Perez is leading an effort to gain more resources at her university for graduate students who are not legal residents. Tuition at the graduate university runs about $22,000 for a standard course load.
"I would say at the undergraduate level there are systems in place that have been a huge success. In graduate school, the resources are very limited," Denisse Rojas, director of Pre-Health Dreamers, a nonprofit immigrant advocacy group focusing on education issues, told the LA Times.
This attempt at Claremont, a private institution in California, is only the latest in an effort to provide resources for illegal immigrants pursuing higher education. Last year, as many as 22 states offered in-state tuition for illegal residents.
While illegal immigrants are now eligible for in-state tuition in many states, out-of-state and international students are still saddled with high tuition rates. A recent video by Campus Reform demonstrated many of their frustrations when students shredded their "Social Security cards" in order to be eligible for in-state tuition.
Perez is currently studying the effects of deportation on college-educated Mexican nationals, but was unable to interview individuals in Mexico because of her fear that she may not be able to return to the United States.
In 2012, Perez published a how-to guide for finding work following college, specifically for undocumented students.
Pre-Health Dreamers and Iliana Perez did not respond to a request for comment from Campus Reform.
Follow the author of this article on Twitter: @LadyLiberty94
Caitlin Grimes
North Carolina Campus Correspondent


text:
Allum is een verdwenen nederzetting op Terschelling waarvan de naam voorkomt in enkele middeleeuwse geschriften. Voor zover bekend moet de nederzetting gelegen hebben ten zuiden van Hee op Terschelling, provincie Friesland (Nederland), in het gebied van de huidige Waddenzee. Vermoed wordt dat het lag op een inmiddels verdwenen buitendijks gebiedje, dat de naam de Oele droeg. In deze omgeving is nog steeds de veldnaam Allumer finne bekend. Allum bestond in 1630 al niet meer.

Middeleeuwse geografie
Geografie van Terschelling
Geschiedenis van Terschelling
Verdwenen plaats in Friesland

redpajama_set_name: RedPajamaWikipedia


For more than 30 years, Eli Broad College of Business' Paulette Stenzel, professor of sustainability and international business law, has given Broad students a grassroots perspective of doing global business by opening the college's doors to some of the world's most innovative and sustainable companies.
Paulette Stenzel has spent much of her life in search of business concepts that align with her personal beliefs. She has traveled across numerous countries strengthening connections between the consumers and the producers. Over the years, she has become invested in finding ways to ensure workers have a say in their labor conditions, and to reinforce transparency in company practices.
Today, Stenzel is an advocate of the fair trade business and a stakeholder within the network of people that operate within these companies.
Through her efforts to promote fair trade Stenzel met Jose Perez Vasquez, a founder of Maya Vinic Coffee Cooperative of Chiapas, Mexico. Maya Vinic comprises more than 500 indigenous coffee farmers who produce more than 400 kilograms of coffee annually.
As part of her global business law course, Stenzel invited Vasquez to speak to Broad students for a real-time example of the fair trade practices they've read about in textbooks. Vasquez shared his story of 47 members of his community caught in political cross-fire and how the aftermath of the tragedy was founding Maya Vinic Coffee Cooperation.
Benefit Corporations, or B-Corporations, like fair trade, focus on the triple bottom line of business because it gives transparency to a company's economic, social and environmental practices. However, it requires legal accountability of companies to continue fair and transparent practices. "It's not always about making shareholders money," said Stenzel, "It's about the human connection and paying fair wages for the workers," she said.
Higher Grounds Trading Company is a certified B-Corporation based in Traverse City, Michigan, where Maya Vinic sells their coffee as fair trade. The co-founder of Higher Grounds, Chris Treter, first met Vasquez when Treter was an intern studying abroad in Chiapas. Higher Grounds places QR codes on product packaging so customers can track the travel and origin to which farmer harvested the coffee. "It is all about the human connection," said Stenzel.
While committed to educating students on sustainable and fair business practices, Stenzel stresses the idea of making sure students find the connection between what companies do and how they themselves use products. "Live your values, be fair to other people. You can practice your values in business," said Stenzel.


text:
Action-from-the-Youth-Americas-Cup---Gilles-Martin-Raget-pic
Red Bull Youth America's Cup final eight decided
Day Two of the Pool A Qualifiers in the Red Bull Youth America's Cup once again belonged to Sweden's Artemis Youth Racing as they topped the tables with two wins.
Heading into the day the Swedish team once again grabbed the headlines with two victories from day two's three races to reaffirm their position at the top of the standings, qualifying comfortably with 55 points.
Joining them in the final three qualification spots from Pool A are SVB Team Germany, Switzerland's Team Tilt and Team France Jeune.
The result today means that Youth Vikings Denmark and Kaijin Team Japan from Pool A join Next Generation USA and Austria's Candidate Sailing Team from Pool B in suffering elimination from the 2017 Red Bull Youth America's Cup.
The results also mean that the identities of the eight teams which will contest the 2017 Red Bull Youth America's Cup are now known. The final field is Artemis Youth Racing, SVB Team Germany, Switzerland's Team Tilt, Team France Jeune, Land Rover BAR Academy, Spanish Impulse by IBEROSTAR, NZL Sailing Team and Bermuda's hometown heroes, TeamBDA.
Red Bull Youth America's Cup Sport Directors Roman Hagara and Hans Peter Steinacher gave their thoughts on the latest day of competition, and what lies ahead in the Finals on 20th and 21st June. Hagara started, saying, "The teams were faced with yet more light winds today, but they dealt with the conditions well and we managed to complete the race program today.
"It's tough work racing in these conditions – it's hot, there's little breeze and you're having to work incredibly hard on the racecourse to maximise every opportunity that comes your way.
"Now we know which teams will be in the Finals next week, and I think we're in for a couple of days of thrilling racing. Hopefully we have some better winds, but they'll go out and deal with whatever is in front of them, and I know that whichever team comes out on top, they'll deserve the spotlight that will be shining upon them."
Hans Peter Steinacher added, "I've been extremely impressed with how all the teams have come into this regatta. There are a few teams emerging as serious contenders, but the margins are so fine, there's really so little to choose between them, I don't think anyone can predict who will come out on top at the end of the series.
"One thing we can predict for sure is how the Finals on Tuesday and Wednesday will captivate Bermuda. TeamBDA did so well to get into the Finals, and now it's all anyone on the island wants to talk about. The guys in that team are doing an incredible job, and now they also have to deal with the hopes and dreams of their entire home country.
"Having seen them in action, both on and off the water, I don't think for a second that will faze them, in fact I think it'll help them as the support is going to be incredible."
Pool A Race 4
Building on their impressive displays on day one of Pool A qualifying, Gustav Petterson's Artemis Youth Racing solidified their position as Pool A leaders with a third overall victory in the qualifying stages to kick off the day in perfect fashion.
It was Max Kohlhoff's SVB Team Germany who had led the way for the majority of the race, but a late final mark surge from the Swedish team ultimately saw them snatch victory ahead of the Germans.
Behind that scrap, a close-fought battle between the four other Pool A teams saw an improved performance from Youth Vikings Denmark help them clinch a morale-boosting third place finish, ahead of Team France Jeune, Switzerland's Team Tilt and Kaijin Team Japan.
Artemis Youth Racing's dominance of Pool A continued in race five, as skipper Gustav Peterson again mastered the light conditions to seal a second successive victory.
It proved a much closer battle for the next three places, with Sebastian Schneiter's Team Tilt securing second place, ahead of SVB Team Germany and Youth Vikings Denmark.
Further back, Ibuki Koizumi's Kaijin Team Japan crossed the finish line in fifth place while somewhat surprisingly, Team France Jeune, who struggled throughout the race, were forced to settle for sixth place.
The results after race five meant that Artemis Youth Racing and SVB Team Germany had secured their places in the Red Bull Youth America's Cup Final. Their progression meant that there were only two spaces left in the Finals for the remaining four teams to fight for ahead of the sixth and final Pool A Qualifiers race.
The final race of the Pool A Qualifiers saw Switzerland's Team Tilt come out on top in a dramatic late tussle with Youth Vikings Denmark, clinching a vital victory which sealed their progression into the Red Bull Youth America's Cup Finals.
Behind them, Team France Jeune also confirmed their qualification with a third place finish, ahead of Artemis Youth Racing and SVB Team German, who had both already sealed their places in the Finals ahead of Race six.
The final race results meant that Artemis Youth Racing, SVB Team Germany, Switzerland's Team Tilt and Team France Jeune have sealed their progression to join NZL Sailing Team, Land Rover BAR Academy, Spanish Impulse by IBEROSTAR and local favourites Team BDA, in the Red Bull Youth America's Cup Final, taking place on June 20-21.
Meanwhile, Youth Vikings Denmark and Kaijin Team Japan joined Next Generation USA and Austria's Candidate Sailing Team, from Pool B, in suffering elimination
Pool A Race Results
Race 4:
1. Artemis Youth Racing (10pts)
2. SVB Team Germany (9pts)
3. Youth Vikings Denmark (8pts)
4. Team France Jeune (7pts)
5. Team Tilt (SUI) (6pts)
6. Kaijin Team Japan (5pts)
1. Artemis Racing (10pts)
1. Team Tilt (SUI) (10pts)
4. Artemis Youth Racing (7pts)
Pool A Final Standings
2. SVB Team Germany (47pts)
3. Team Tilt (SUI) (46 pts)
4. Team France Jeune (46pts)
5. Youth Vikings Denmark (43pts)
6. Kaijin Team Japan (33pts)
Red Bull Youth America's Cup media hub: www.redbullsailingnewsroom.com
35th America's Cup
End of an era as McConaghy's Jono Morris heads to MalaysiaClose racing in Hobart's SB20 mid-winters


Travel adjustments are commonplace on today's suspension forks: If there isn't an external adjustment, you can usually rearrange, add or remove spacers inside the fork to change travel. But these days, travel adjustments might not be only about travel: With the increase in popularity of 29-inch-wheel bikes, limiting handlebar height is becoming a greater challenge on XC and trail bikes—especially for smaller riders. So, does it make sense to use a fork's travel adjustments to improve bike fit? We went to Brandon Sloan, Specialized's mountain bike product manager, for answers.
Mountain Bike: In what circumstances will you use fork height and travel to affect bike fit?
Brandon Sloan: We use fork travel (or lack of it) in some cases to help the bike fit certain riders better. Some of our RockHopper and HardRock models have different fork travels for different size bikes (smaller hardtail frames run 80mm, while medium and large go to 100mm). This is to help keep standover as low as possible and the bar low for smaller riders.
MB: Do you do anything special or have any limitations regarding fork travel and fit when designing 29ers?
BS: On 29ers, fork travel is a bit tricky to deal with, as the bar can get too high really quick. For the past two years, our high-end 29er has had a bit less travel to help keep the front end low (90mm on 26- inch hardtails, 80mm on 29-inch hardtails). That short-changed the 29ers, so we're trying to level it out in 2010. The problem is that you have to make shorter head tubes to keep the bar low. A short head tube does not help the frame's strength and stiff ness (it is actually the opposite), so there are limits to how short a head tube we can pull off without overbuilding the rest of the frame to make up for it.
MB: For riders who decrease fork travel to get a lower handlebar height, what effect does that have on a bike's geometry?
BS: If you adjust your fork down from the original spec the geometry was based around, you will get a steeper head angle (20mm lower will increase the angle by about 1 degree...lots for a hardtail); a lower bottom bracket (how much lower depends on frame size); a steeper seat tube angle (so the rider might need to readjust seat position); and a lower handlebar height.
MB: What effect does that have on handling?
BS: A steeper head angle makes the bike's steering quicker and less stable at speed. The lower bottom bracket reduces ground clearance but makes for a nice-turning bike. The steeper seat-tube angle moves the saddle farther forward over the BB, and the lower bar height can make the riding position a little more aggressive. Sometimes weight can be too far forward, but if people are making a travel change [to adjust fit], the bar was already too high.
MB: How much can a rider change fork travel before the effects on handling and fit become too great?
BS: I think geometry is the most important part of a bike. The nicest-riding fork or rear suspension will be overshadowed by jacked-up geometry. Changing travel to adjust geometry can be done as long as the rider isn't happy with the original geometry. Ten millimeters might be okay for some riders, but 20mm is way too much (again, unless they don't like the stock geometry).


text:
\section{Introduction}
\seclab{intro}
Inverse problems are ubiquitous in science and engineering. Perhaps
the most popular family of inverse problems is to determine a set of
parameters (or a function) given a set of indirect observations, which
are in turn provided by a parameter-to-observable map plus observation
uncertainties. For example, if one considers the problem of
determining the heat conductivity of a thermal fin given
measured temperature at a few locations on the thermal fin, then: i)
the desired unknown parameter is the distributed heat conductivity,
ii) the observations are the measured temperatures, iii) the
parameter-to-observable map is the mathematical model that describes
the temperature on the thermal fin as a function of the heat
conductivity; indeed the temperature distribution is a solution of an
elliptic partial differential equation (PDE) whose coefficient is the
heat conductivity, and iv) the observation uncertainty is due to the
imperfection of the measurement device and/or model
inadequacy.

The Bayesian inversion framework refers to a mathematical method that
allows one to solve statistical inverse problems taking into account
all uncertainties in a systematic and coherent manner. The Bayesian
approach does this by reformulating the inverse problem as a problem
in statistical inference, incorporating uncertainties in the
observations, the parameter-to-observable map, and prior information
on the parameter. In particular, we seek a statistical description of
all possible (set of) parameters that conform to the available prior knowledge
and at the same time are consistent with the observations. The
solution of the Bayesian framework is the so-called posterior measure
that encodes the degree of confidence on each set of parameters as the
solution to the inverse problem under consideration.

Mathematically the posterior is a surface in high dimensional
parameter space. The task at hand is therefore to explore the
posterior by, for example, characterizing the mean, the covariance,
and/or higher moments. The nature of this task is to compute high
dimensional integrals for which most contemporary methods are
intractable. Perhaps the most general method to attack these problems
is the Markov chain Monte Carlo (MCMC) method which shall be introduced in subsequent sections.

Let us now summarize the content of the paper. We start with the
description of the statistical inverse problem under consideration in
Section \secref{infiniteBayes}. It is an inverse steady state heat
conduction governed by elliptic PDEs. We postulate a Gaussian measure
prior on the parameter space to ensure that the inverse problem is
well-defined. The prior itself is a well-defined object whose
covariance operator is the inverse of an elliptic differential operator and
with the mean function living in the Cameron-Martin space of the
covariance. The posterior is given by its Radon-Nikodym derivative
with respect to the prior measure, which is proportional to the
likelihood. Since the RMHMC simulation method requires the gradient, Hessian, and the
derivative of the Fisher information operator, we discuss, in some
depth, how to compute the derivatives of the potential function
(the misfit functional) with PDE constraints efficiently using the adjoint
technique in Section \secref{adjoint}. In particular, we define a
Fisher information operator and show that it coincides with the
well-known Gauss-Newton Hessian of the misfit. We next present a
discretization scheme for the infinite Bayesian inverse problem in
Section \secref{FEM}. Specifically, we employ a standard continuous
$H^1$-conforming finite element (FEM) method to discretize both the
likelihood and the Gaussian prior. We choose to numerically compute
the truncated Karhunen-Lo\`eve expansion which requires one to solve an
eigenvalue problem with fractional Laplacian. In order to accomplish
this task, we use a matrix transfer technique (MTT) which leads to a
natural discretization of the Gaussian prior measure.  In Section
\secref{HMC}, we describe the Riemannian manifold Hamiltonian Monte
Carlo (RMHMC) and its variants at length, and its application to our
Bayesian inverse problem. Section \secref{lowrank} presents a low rank
approach to approximate the Fisher information matrix and its inverse
efficiently. This is possibly due to the fact that the Gauss-Newton
Hessian, and hence the Fisher information operator, is a compact
operator. Various numerical results supporting our proposed approach
are presented in Section \secref{results}. We begin this section with
an extensive study and comparison of Riemannian manifold MCMC
methods for problems with two parameters, and end the section with
$1025$-parameter problem. Finally, we conclude the paper in Section
\secref{conclusions} with a discussion on future work.

\section{Problem statement}
\seclab{infiniteBayes} In order to clearly illustrate the challenges arising in PDE-constrained inverse problems for MCMC based Bayesian inference, we consider
the following heat conduction problem governed by an elliptic partial
differential equation in the open and bounded domain $\Omega \subset \R^n$:
\begin{align*}
-\Div\LRp{e^\u\Grad \w} &= 0 & \text{ in } \Omega \\
-e^\u\Grad \w \cdot \mb{n} &= Bi \,u &\text{ on } \partial \Omega \setminus \Gamma_{R}, \\
 -e^\u\Grad \w \cdot \mb{n} &= -1 &\text{ on } \Gamma_{R},
 \end{align*}
where $\w$ is the forward state, $\u$ the logarithm of distributed thermal conductivity on $\Omega$, $\mb{n}$ the unit outward normal on
$\pOmega$, and $Bi$ the Biot number.

In the forward problem, the task is to solve for the temperature
distribution $\w$ given a description of distributed parameter $\u$.
In the inverse problem, the task is to reconstruct $\u$ given some
available observations, e.g, temperature observed at some
parts/locations of the domain $\Omega$. We initially choose to cast
the inverse problem in the framework of PDE-constrained
optimization. To begin, let us consider the following additive
noise-corrupted pointwise observation model\footnote{We assume the
  forward state $\w$ is sufficiently regular, i.e. $\w \in H^{s}, s >
  n/2$, so that $w$ is, by the virtue of the Sobolev embedding
  theorem, continuous, and therefore it is meaningful to measure $w$
  pointwise.}
\begin{equation}\eqnlab{pointwiseObs}
\d_j := \w\LRp{\x_j} + \eta_j, \quad j = 1,\hdots,K,
\end{equation}
where $K$ is the total number of observation locations,
$\LRc{\mb{x}_j}_{j=1}^K$ the set of points at which $\w$ is observed,
$\eta_j$ the additive noise, and $\d_j$ the actual noise-corrupted
observations. In this paper we work with synthetic observations and
hence there is no model inadequacy in
\eqnref{pointwiseObs}. Concatenating all the observations, one can
rewrite \eqnref{pointwiseObs} as
\begin{equation}
\eqnlab{observation}
\db := \mc{G}\LRp{\u} + \etab,
\end{equation}
with $\mc{G} := \LRs{\w\LRp{\mb{x}_1}, \hdots,\w\LRp{\mb{x}_K}}^T$
denoting the map from the distributed parameter $\u$ to the
noise-free observables, $\etab$ being random numbers normally distributed by
$\GM{0}{\L}$ with bounded covariance matrix $\L$, and $\db =
\LRs{\d_1,\hdots,\d_K}^T$. For simplicity, we take $\L = \sigma^2\I$,
where $\I$ is the identity matrix.

Our inverse problem can be now formulated as
\begin{align}
&\min_{\u}  \J\LRp{\u,\db} := \frac{1}{2}\snor{\db-\mc{G}\LRp{\u}}_\L^2 =
\frac{1}{2\sigma^2} \sum_{j = 1}^K\LRp{\w\LRp{\x_j}-\d_j}^2
\eqnlab{Cost}
\end{align}
\SubjectTo
\begin{subequations}
\eqnlab{forwardEqn}
\begin{align}
-\Div\LRp{e^\u\Grad \w} &= 0 & \text{ in } \Omega, \eqnlab{forward} \\
-e^\u\Grad \w \cdot \mb{n} &= Bi \,u &\text{ on } \partial \Omega \setminus \Gamma_{R}, \eqnlab{forwardbcRobin}\\
 -e^\u\Grad \w \cdot \mb{n} &= -1 &\text{ on } \Gamma_{R},\eqnlab{forwardbcNeumann}.
\end{align}
\end{subequations}
where $\snor{\cdot}_\L:= \snor{\L^{-\half}\cdot}$ denotes the weighted
Euclidean norm induced by the canonical inner product
$\LRp{\cdot,\cdot}$ in $\R^K$. This optimization problem is however
ill-posed. An intuitive reason is that the dimension of observations
$\db$ is much smaller than that of the parameter $\u$, and hence they
provide limited information about the distributed parameter $\u$. As a
result, the null space of the Jacobian of the parameter-to-observation map
$\mc{G}$ is non-empty. Indeed, we have shown that the Gauss-Newton
approximation of the Hessian (which is the square of this Jacobian,
and is also equal to the full Hessian of the data misfit $\mathcal{J}$
evaluated at the optimal parameter) is a compact operator
\cite{Bui-ThanhGhattas12a, Bui-ThanhGhattas12, Bui-ThanhGhattas12f},
and hence its range space is effectively finite-dimensional.

One way to overcome the ill-posedness is to use {\em Tikhonov
  regularization} (see, e.g., \cite{Vogel02}), which proposes to
augment the cost functional \eqnref{Cost} with a quadratic term, i.e.,
\begin{equation}
\tilde{\mc{J}} := \frac{1}{2}\snor{\db-\mc{G}\LRp{\u}}_\L^2 +
\frac{\kappa}{2}\nor{R^{1/2}\u}^2,
\eqnlab{CostRegularized}
\end{equation}
where $\kappa$ is a regularization parameter, $R$ some regularization
operator, and $\nor{\cdot}$ some appropriate norm.  This method is a
representative of deterministic inverse solution techniques that
typically do not take into account the randomness due to measurements
and other sources, though one can equip the deterministic solution
with a confidence region by post-processing (see, e.g.,
\cite{Vexler04} and references therein). It should be pointed out that
if the regularization term is replaced by the Cameron-Martin norm of
$\u$ (the second term in \eqnref{MAP}), the Tikhonov solution is in
fact identical to the maximum a posteriori point in
\eqnref{MAP}. However, such a point estimate is insufficient for the
purpose of fully taking the randomness into account.

In this paper, we choose to tackle the ill-posedness using a {\em
  Bayesian} framework \cite{Franklin70,
  LehtinenPaivarintaSomersalo89, Lasanen02, Stuart10, Piiroinen05}. We
seek a statistical description of all possible $\u$ that conform to
some prior knowledge and at the same time are consistent with the
observations. The Bayesian approach does this by reformulating the
inverse problem as a problem in {\em statistical inference},
incorporating uncertainties in the observations, the forward map
$\mc{G}$, and prior information. This approach is appealing since it
can incorporate most, if not all, kinds of randomness in a systematic
manner.  To begin, we postulate a Gaussian measure $\mu :=
\GM{\u_0}{\alpha^{-1}\mc{C}}$ on $\u$ in $\Ltwo$ where
\[
\mc{C} := \LRp{I - \Delta}^{-s} =: \mc{A}^{-s}
\]
with the domain of definition
\[
D\LRp{\mc{A}} := \LRc{u \in H^2\LRp{\Omega}: \pp{u}{\mb{n}} = 0 \text{ on } \pOmega},
\]
where $H^2\LRp{\Omega}$ is the usual Sobolev space. Assume that the
mean function $\u_0$ lives in the Cameron-Martin space of $\mc{C}$,
then one can show (see \cite{Stuart10}) that the measure $\mu$ is
well-defined when $s > n/2$ ($d$ is the spatial dimension), and in
that case, any realization from the prior distribution $\mu$ is almost
surely in the H\"older space $\X := C^{0,\beta}\LRp{\Omega}$ with $0 <
\beta < s/2$. That is, $\mu\LRp{X} = 1$, and the
Bayesian posterior measure $\nu$ satisfies the Radon-Nikodym
derivative
\begin{equation}
\eqnlab{RadonNikodym}
\pp{\nu}{\mu}\LRp{\u|\db} \sim \exp\LRp{-\mc{J}\LRp{\u,\db}} = \exp\LRp{-\half\snor{\db -
    \mc{G}\LRp{\u}}^2_\L},
\end{equation}
if $\mc{G}$ is a continuous map from $\X$ to $\R^K$. Note that the
Radon-Nikodym derivative is proportional to the the likelihood defined by
\[
\like \sim \exp\LRp{-\mc{J}\LRp{\u,\db}}.
\]

The maximum a posteriori (MAP) point is
defined as
\begin{equation}
\eqnlab{MAP}
\u^{MAP} := \arg\min_\u \mc{J}\LRp{\u, \db} :=\half\snor{\db -
    \mc{G}\LRp{\u}}^2_\L + \frac{\alpha}{2} \nor{\u}^2_{\mc{C}},
\end{equation}
where $\nor{\cdot}_{\mc{C}} :=
\nor{\mc{C}^{-\half}\cdot}$ denotes the weighted
$\Ltwo$ norm induced by the $\Ltwo$ inner product $\LRa{\cdot,\cdot}$.

\section{Adjoint computation of gradient, Hessian, and the third derivative tensor}
\seclab{adjoint} In this section, we briefly present the adjoint
method to efficiently compute the gradient, Hessian,
and the third derivative of the cost functional \eqnref{Cost}. We
start by considering the weak form of the (first order) forward
equation \eqnref{forwardEqn}:
\begin{equation}
\eqnlab{Forwardw}
\iOm{e^\u\Grad\w \cdot \Grad\lamh} + \iGb{Bi\, \w \lamh} = \iGs{\lamh},
\end{equation}
with $\lamh$ as the test function. Using the standard reduced space
approach (see, e.g., a general discussion in \cite{NocedalWright06}
and a detailed derivation in \cite{Bui-ThanhGhattas13}) one can show
that the gradient $\Grad\mc{J}\LRp{\u}$, namely the Fr\'echet
derivative of the cost functional $\mc{J}$, acting in any direction
$\uo$ is given by
\begin{equation}
\eqnlab{Gradient}
\LRa{\Grad\mc{J}\LRp{\u},\uo} =
\iOm{\uo e^{\u}\Grad\w\cdot\Grad\lambda},
\end{equation}
where the (first order) adjoint state $\lambda$ satisfies the adjoint equation
\begin{equation}
\eqnlab{Adjointw}
\iOm{e^\u\Grad\lambda \cdot \Grad\wh} + \iGb{Bi\, \lambda \wh} =
-\frac{1}{\sigma^2}\sum_{j=1}^K\LRp{\w\LRp{\x_j}-\d_j}\wh\LRp{\x_j},
\end{equation}
with $\wh$ as the test function.  On the other hand, the Hessian, the
Fr\'echet derivative of the gradient, acting in directions $\uo$ and
$\ut$ (superscript ``2'' means the second variation direction) reads
\begin{equation}
\eqnlab{FullHessian}
\LRa{\LRa{\Grad^2\mc{J}\LRp{\u},\uo},\ut} =
\iOm{\uo e^\u \Grad\w\cdot\Grad\lamt} + \iOm{\uo e^\u\Grad\wt\cdot\Grad\lambda} +
\iOm{\uo\ut e^\u\Grad\w\cdot\Grad\lambda},
\end{equation}
where the second order forward state $\wt$ obeys the second order forward equation
\begin{equation}
\eqnlab{Iforward}
\iOm{e^\u\Grad\wt \cdot \Grad\lamh} + \iGb{Bi\, \wt \lamh} = - \iOm{\ut e^\u\Grad\w \cdot \Grad\lamh},
\end{equation}
and the second order  adjoint state $\lamt$ is governed by the second order adjoint equation
\begin{equation}
\eqnlab{Iadjoint}
\iOm{e^\u\Grad\lamt \cdot \Grad\wh} + \iGb{Bi\, \lamt \wh} =
-\frac{1}{\sigma^2}\sum_{j=1}^K\wt\LRp{\x_j}\wh\LRp{\x_j} - \iOm{\ut e^\u\Grad\lambda \cdot \Grad\wh}.
\end{equation}

We define the generalized Fisher information operator\footnote{Note that the
  Fisher information operator is typically defined for finite dimensional settings in which it is a matrix.}  acting in directions $\uo$ and $\ut$ as
\begin{equation}
\eqnlab{Fisher}
\LRa{\LRa{G\LRp{\u},\uo },\ut} := \Ex_\like\LRs{\LRa{\LRa{\Grad^2\mc{J}\LRp{\u},\uo},\ut}},
\end{equation}
where the expectation is taken with respect to the likelihood---the
distribution of the observation $\db$. Now, substituting
\eqnref{FullHessian} into \eqnref{Fisher} and assuming that the
integrals/derivatives can be interchanged we obtain
\begin{align*}
\LRa{\LRa{G\LRp{\u},\uo },\ut} &= \iOm{\uo e^\u \Grad\w\cdot\Grad\Ex_\like\LRs{\lamt}} + \iOm{\uo e^\u\Grad\wt\cdot\Grad\Ex_\like\LRs{\lambda}} \\
&+
\iOm{\uo\ut e^\u\Grad\w\cdot\Grad\Ex_\like\LRs{\lambda}},
\end{align*}
where we have used the assumption that the parameter $\u$ is
independent of observation $\db$ and the fact that $\w$ and $\wt$ do
not depend on $\db$. The next step is to compute
$\Grad\Ex_\like\LRs{\lambda}$ and $\Ex_\like\LRs{\lamt}$. To begin,
let us take the expectation the first order adjoint equation \eqnref{Adjointw} with respect to $\like$ to arrive at
\begin{multline*}
\iOm{e^\u\Grad\Ex_\like\LRs{\lambda} \cdot \Grad\wh} + \iGb{Bi\, \Ex_\like\LRs{\lambda} \wh} =
\\-\frac{1}{\sigma^2}\sum_{j=1}^K\Ex_\like\LRs{\w\LRp{\x_j}-\d_j}\wh\LRp{\x_j} = 0,
\end{multline*}
where the second equality is obtained from \eqnref{pointwiseObs} and
the assumption $\eta_j \sim \GM{0}{\sigma^2}$. We conclude that
\begin{equation}
\eqnlab{ExAdjoint}
\Grad\Ex_\like\LRs{\lambda}  = 0.
\end{equation}
On the other hand, if we take the expectation of the second order adjoint equation
\eqnref{Iadjoint} and use \eqnref{ExAdjoint} we have
\begin{equation}
\eqnlab{ExIadjoint}
\iOm{e^\u\Grad\Ex_\like\LRs{\lamt} \cdot \Grad \wh} + \iGb{Bi\, \Ex_\like\LRs{\lamt} \wh} =
-\frac{1}{\sigma^2}\sum_{j=1}^K\wt\LRp{\x_j}\wh\LRp{\x_j}.
\end{equation}
Let us define
\[
\lamtt := \Ex_\like\LRs{\lamt},
\]
then \eqnref{ExIadjoint} becomes
\begin{equation}
\eqnlab{ExIadjointF}
\iOm{e^\u\Grad\lamtt \cdot \Grad\wh} + \iGb{Bi\, \lamtt \wh} =
-\frac{1}{\sigma^2}\sum_{j=1}^K\wt\LRp{\x_j} \wh\LRp{\x_j}.
\end{equation}

As a result, the Fisher information operator acting along directions
$\uo$ and $\ut$ reads
\begin{equation}
\eqnlab{FisherF}
\LRa{\LRa{G\LRp{\u},\uo },\ut} = \iOm{\uo e^\u \Grad\w\cdot\Grad\lamtt},
\end{equation}
where $\lamtt$ is the solution of \eqnref{ExIadjointF}, a variant of
the second order adjoint equation \eqnref{Iadjoint}. The Fisher
information operator therefore coincides with the Gauss-Newton Hessian of the cost functional \eqnref{Cost}.

The procedure for computing the gradient acting on an arbitrary
direction is clear. One first solves the first order forward equation
\eqnref{Forwardw} for $\w$, then the first order adjoint
\eqnref{Adjointw} for $\lambda$, and finally evaluate
\eqnref{Gradient}. Similarly, one can compute the Hessian (or the
Fisher information operators) acting on two arbitrary directions by
first solving the second order forward equation \eqnref{Iforward} for
$\wt$, then the second order adjoint equation \eqnref{Iadjoint} (or
its variant \eqnref{ExIadjointF}) for $\lamt$ (or $\lamtt$), and
finally evaluating \eqnref{FullHessian} (or \eqnref{FisherF}).

One of the main goals of the paper is to study the Riemann manifold
Hamiltonian Monte Carlo method in the context of Bayesian inverse problems
governed by PDEs. It is therefore essential to compute the derivative
of the Fisher information operator. This task is obvious for problems
with available closed form expressions of the likelihood and the prior,
but it is not so for those governed by PDEs. Nevertheless, using the
adjoint technique we can compute the third order derivative tensor
acting on three arbitrary directions with three extra PDE solves, as
we now show. To that end, recall that the Fisher information operator
acting on directions $\uo$ and $\ut$ is given by \eqnref{FisherF}. The
Fr\'echet derivative of the Fisher information operator along the additional
direction $\uth$ (superscript ``3'' means the third variation direction) is given by
\begin{align}
&\LRa{\LRa{\LRa{\T\LRp{\u}, \uo},\ut}, \uth} := \LRa{\Grad\LRa{\LRa{G\LRp{\u}, \uo},\ut}, \uth} \nonumber \\&= \iOm{\uo\uth e^\u\Grad\w \cdot \Grad\lamtt} + \iOm{\uo e^\u \Grad \wth \cdot \Grad\lamtt} + \iOm{\uo e^\u \Grad\w\cdot\Grad\lamtth}
\eqnlab{tensor},
\end{align}
where $\wth$, $\lamtth$ are the variation of $\w$ and $\lamtt$ in the
direction $\uth$, respectively. One can show that $\uth$ satisfies
another second order forward
equation
\begin{equation}
\eqnlab{Iforwardn}
\iOm{e^\u\Grad\wth \cdot \Grad \lamh} + \iGb{Bi\, \wth \lamh} = - \iOm{\uth e^\u\Grad\w \cdot \Grad\lamh}.
\end{equation}
Similarly, $\lamtth$ is the solution of the third order adjoint equation
\begin{equation}
\eqnlab{IIadjoint}
\iOm{e^\u\Grad\lamtth \cdot \Grad \wh} + \iGb{Bi\, \lamtth \wh} =
-\frac{1}{\sigma^2}\sum_{j=1}^K\wtth\LRp{\x_j} \wh\LRp{\x_j}
 - \iOm{\uth e^\u\Grad\lamtt \cdot \Grad\wh},
\end{equation}
and $\wtth$, the variation of $\ut$ in direction $\uth$, satisfies the
following third order forward equation
\begin{align}
&\iOm{e^\u\Grad\wtth \cdot \Grad \lamh} + \iGb{Bi\, \wtth \lamh} = \nonumber\\
&-\iOm{\uth e^\u\Grad\wt \cdot \Grad \lamh}
- \iOm{\uth\ut e^\u\Grad\w \cdot \Grad\lamh}
- \iOm{\ut e^\u\Grad\wth \cdot \Grad\lamh}.
\eqnlab{IIforward}
\end{align}
Note that it would have required four extra PDE solves if one computes
the third derivative of the full Hessian \eqnref{FullHessian}.

It is important to point out that the operator $\T$ is only symmetric
with respect to $\uo$ and $\ut$ since the Fisher information is
symmetric, but not with respect to $\uo$ and $\uth$ or $\ut$ and
$\uth$. The full symmetry only holds for the derivative of the full
Hessian, that is, the true third derivative of the cost functional.

\section{Discretization}
\seclab{FEM} As presented in Section \secref{infiniteBayes}, we view
our inverse problem from an infinite dimensional point of view. As
such, to implement our approach on computers, we need to discretize
the prior, the likelihood and hence the posterior. We choose to use the finite element method.
In particular, we employ the
standard $H^1\LRp{\Omega}$ finite element method (FEM) to discretize
the forwards and adjoints (the
likelihood), and the operator $\mc{A}$ (the prior). It should be
pointed out that the Cameron-Martin space can be shown (see, e.g.,
\cite{Stuart10}) to be a subspace of the usual fractional Sobolev
space $H^s\LRp{\Omega}$, which is in turn a subspace of
$H^1\LRp{\Omega}$. Thus, we are using a non-conforming FEM approach
(outer approximation). For convenience, we further assume that
the discretized state and parameter live on the same finite
element mesh. Since FEM approximation of elliptic operators is standard (see, e.g., \cite{Ciarlet78}), we will not discuss it here. Instead, we describe the matrix transfer technique (see, e.g, \cite{IlicLiuTurnerEtAl05} and the references therein) to
discretize the prior.

Define $\Q := \mc{C}^{1/2} = \mc{A}^{-s/2}$, then the eigenpairs
$\LRp{\lambda_i,\v_i}$ of $\Q$ define the Karhunen-Lo\`eve (KL)
  expansion of the prior distribution as
\[
\u = \u_0+ \frac{1}{\sqrt{\alpha}}\sum_{i=1}^\infty a_i \lambda_i \v_i,
\]
where $a_i \sim \GM{0}{1}$. We need to solve
\[
\Q \v_i = \lambda_i\v_i,
\]
or equivalently
\begin{equation}
\eqnlab{eigenProblem}
\mc{A}^{s/2} \v_i = \frac{1}{\lambda_i}\v_i.
\end{equation}

To solve \eqnref{eigenProblem} using the matrix transfer technique (MTT),
let us denote by $\M$ the mass matrix, and $\K$ the stiffness matrix
resulting from the discretization of the Laplacian $\Delta$. The
representation of $\mc{A}$ in the finite element space (see, e.g., \cite{Simpson08} and the references therein)  is given by
\[
\A := \M^{-1}\K + \I.
\]
Let bold symbols denote the corresponding vector of FEM nodal values,
e.g., $\ub$ is the vector containing all FEM nodal values of
$\u$. If we define $\LRp{\sigma_i,\vb_i}$ as
eigenpairs for $\A$, i.e,
\begin{equation}
\eqnlab{eigenmatrix}
\A \vb_i = \sigma_i \vb_i, \text{ or } \A \Vb = \boldsymbol{\Sigma}\Vb
\end{equation}
where $\vb_i^T\M\vb_j = \delta_{ij}$, and hence $\Vb^{-1} = \Vb^T\M$,
$\delta_{ij}$ is the Kronecker delta function, and $\boldsymbol{\Sigma}$ is the
diagonal matrix with entries $\sigma_i$. Since $\A$ is similar to $\M^{-\half}\LRp{\K +
  \M}\M^{-\half}$, a symmetric positive definite matrix, $\A$ has
positive eigenvalues.  Using MTT method, the matrix representation of
\eqnref{eigenProblem} reads
\[
\A^{s/2}\vb_i = \frac{1}{\lambda_i}\vb_i,
\]
where
\[
\A^{s/2} := \Vb\boldsymbol{\Sigma}^{s/2}\Vb^{-1}.
\]
It follows that
\[
\lambda_i = \sigma^{-s/2}_i.
\]
The Galerkin FEM approximation of the prior via truncated KL expansion reads
\begin{equation}
\eqnlab{uKLtruncated}
\ub = \ub_0+ \frac{1}{\sqrt{\alpha}}\sum_{i=1}^{N}a_i\lambda_i\vb_i,
\end{equation}
with $\ub$ as the FEM nodal value of the approximate prior sample $\u$
and $N$ as the number of FEM nodal points. Note that for ease in
writing, we have used the same notation $\u$ for both infinite
dimensional prior sample and its FEM approximation. Since $\u
\in \Ltwo$, $\ub$ naturally lives in $\R^N_\M$, the Euclidean
space with weighted inner product $\LRp{\cdot,\cdot}_\M := \LRp{\cdot,\M\cdot}$.

A question arises: what is the distribution of $\ub$? Clearly
$\ub$ is a Gaussian with mean $\ub_0$ since $a_i$ are. The
covariance matrix $\C$ for $\ub$ is defined by
\[
\LRp{\zb, \C\yb}_\M := \Ex\LRs{\LRp{\ub-\ub_0,\M\zb}\LRp{\ub-\ub_0,\M\yb}} = \frac{1}{\alpha}\zb^T \M \Vb \boldsymbol{\Lambda}^2 \Vb^T\M\yb,
\]
where we have used \eqnref{uKLtruncated} to obtain the second equality and  $\boldsymbol{\Lambda}$ is the diagonal matrix with entries $\Lambda_{ii}
= \lambda_i^{-1}$. It follows that
\begin{equation}
\eqnlab{eqnC}
\C = \frac{1}{\alpha} \Vb \boldsymbol{\Lambda}^2 \Vb^T\M
\end{equation}
as a map from $\R^N_\M$ to $\R^N_\M$, and its inverse can be shown to be
\[
\C^{-1} = \alpha \Vb \boldsymbol{\Lambda}^{-2} \Vb^T\M,
\]
whence the distribution of $\ub$ is
\begin{equation}
\eqnlab{uDist}
\ub \sim \GM{\ub_0}{\alpha \Vb \boldsymbol{\Lambda}^{-2} \Vb^T\M} \sim
\exp\LRs{-\frac{\alpha}{2} \LRp{\ub-\ub_0}^T\M\Vb \boldsymbol{\Lambda}^{-2} \Vb^T\M\LRp{\ub-\ub_0}}.
\end{equation}
As a result, the FEM discretization of the prior can be written as
\begin{align*}
\frac{\alpha}{2}\nor{\u-\u_0}^2_{\mc{C}} &:= \frac{\alpha}{2} \nor{\mc{A}^{s/2}\LRp{\u-\u_0}}^2 \stackrel{\text{MTT}}{\approx}
 \frac{\alpha}{2}\LRp{\ub-\ub_0}^T\M\Vb\boldsymbol{\Lambda}^{-2}\Vb^T\M\LRp{\ub-\ub_0}.
\end{align*}
Thus, the FEM approximation of the posterior is given by
\[
\post \sim \exp\LRp{-\half\snor{\db -
    \mc{G}\LRp{\u}}^2_\L} \times \exp\LRp{-\frac{\alpha}{2}\LRp{\ub-\ub_0}^T\M\Vb\boldsymbol{\Lambda}^{-2}\Vb^T\M\LRp{\ub-\ub_0}}.
\]
The detailed derivation of the FEM approximation of infinite Bayesian
inverse problems in general and the prior in particular will be
presented elsewhere \cite{Bui-Thanh14}.

\section{Riemannian manifold Langevin and Hamiltonian Monte Carlo methods}
\seclab{HMC} In this section we give a brief overview of the MCMC
algorithms that we consider in this work. Some familiarity with the
concepts of MCMC is required by the reader since an introduction to
the subject is out of the scope of this paper.

\subsection{Metropolis-Hastings}
For a random vector $\ub \in \mathbb{R}^N$ with
density $\pi(\ub)$ the Metropolis-Hastings algorithm
employs a proposal mechanism
$q(\ub^{*}|\ub^{t-1})$ and proposed
moves are accepted with probability $\min
\left\{1,\pi(\ub^{*})
q(\ub^{t-1}|\ub^{*})/\pi(\ub^{t-1})q(\ub^{*}|\ub^{t-1})\right\}$. Tuning
the Metropolis-Hastings algorithm involves selecting an appropriate
proposal mechanism. A common choice is to use a Gaussian proposal of
the form $q(\ub^{*}|\ub^{t-1}) =
\mathcal{N}(\ub^{*}|\ub^{t-1},\boldsymbol{\Sigma})$,
where $\mathcal{N}(\cdot|\boldsymbol{\mu},\boldsymbol{\Sigma})$
denotes the multivariate normal density with mean $\boldsymbol{\mu}$
and covariance matrix $\boldsymbol{\Sigma}$.

Selecting the covariance matrix however, is far from trivial in most
of cases since knowledge about the target density is
required. Therefore a more simplified proposal mechanism is often
considered where the covariance matrix is replaced with a diagonal
matrix such as $\boldsymbol{\Sigma}=\epsilon\I$ where the
value of the scale parameter $\epsilon$ has to be tuned in order to
achieve fast convergence and good mixing. Small values of $\epsilon$
imply small transitions and result in high acceptance rates while the
mixing of the Markov Chain is poor. Large values on the other hand,
allow for large transitions but they result in most of the samples
being rejected. Tuning the scale parameter becomes even more difficult
in problems where the standard deviations of the marginal posteriors
differ substantially, since different scales are required for each
dimension, and when correlations between different variables exist. In
the case of PDE-constrained inverse problems in very high dimensions
with strong nonlinear interactions inducing complex non-convex
structures in the target posterior this tuning procedure is typically
doomed to failure of convergence and mixing.

There have been many subsequent developments of this basic algorithm
however the most important with regard to inverse problems, arguably,
is the formal definition of Metropolis Hastings in an infinite
dimensional functional space. One of the main failings of Metropolis
Hastings is the drop-off in acceptance probability as the dimension of
the problem increases. By defining the Metropolis acceptance
probability in the appropriate Hilbert space the acceptance
probability should then be invariant to the dimension of the problem
and this is indeed the case as is described in a number of scenarios
by \cite{CotterRobertsStuartEtAl13}. Furthermore the definition of a Markov chain
transition kernel directly in the Hilbert space which exploits
Hamiltonian dynamics in the proposal mechanism followed in
\cite{OttobrePillaiPinskiEtAl14}.

These are important methodological advances for MCMC applied to
Inverse Problems. As the infinite dimensional nature of the problem is
a fundamental aspect of the problem it is sensible that this
characteristic is embedded in the MCMC scheme. In a similar vein by
noting that the statistical model associated with the specific inverse
problem is generated from an underlying partial differential equation
or system of ordinary differential equations a natural geometric
structure structure on the space of probability distributions is
induced. This structure provides a rich source of model specific
information that can be exploited in devising MCMC schemes that are
informed by the underlying structure of the model itself.

In \cite{GirolamiCalderhead11} a way around this situation was
provided by accepting that the statistical model can itself be
considered as an object with an underlying geometric structure that
could be embedded into the proposal mechanism. A class of MCMC methods
were developed based on the differential geometric concepts underlying
Riemannian manifolds.

\subsection{Riemann Manifold Metropolis Adjusted Langevin Algorithm}
Denoting the log of the target density as $\mathcal{L}(\ub) = \log \pi(\ub)$, the manifold Metropolis Adjusted Langevin Algorithm (mMALA) method, \cite{GirolamiCalderhead11}, defines a Langevin diffusion with stationary distribution $\pi(\ub)$ on the Riemann manifold of density functions with metric tensor $\G(\ub)$. By employing a first order Euler integrator for discretising the stochastic differential equation a proposal mechanism with density
$q(\ub^*|\ub^{t-1}) = \mathcal{N}(\ub^*| \boldsymbol{\mu}(\ub^{t-1},\epsilon),\epsilon^2\G^{-1}(\ub^{t-1}))$
is defined, where $\epsilon$ is the integration step size, a parameter which needs to be tuned, and the $k$th component of the mean function $\boldsymbol{\mu}(\ub,\epsilon)_k$ is
\begin{eqnarray}
\boldsymbol{\mu}(\ub,\epsilon)_k & = & \ub_k +
\frac{\epsilon^2}{2}\left(\G^{-1}(\ub)\nabla_{\ub}\mathcal{L}(\ub)\right)_k
- \epsilon^2
\sum_{i=1}^N\sum_{j=1}^N\G(\ub)_{i,j}^{-1}\Gamma_{i,j}^k \label{eq:meanmMALA}\end{eqnarray}
where $\Gamma_{i,j}^k$ are the Christoffel symbols of the metric in
local coordinates. Note that we have used the Christoffel symbols to
express the derivatives of the metric tensor, and they are computed using the adjoint method presented in Section \secref{adjoint}.

Due to the discretisation error introduced by the first order approximation convergence to the stationary distribution is not guaranteed anymore and thus the Metropolis-Hastings ratio is employed to correct for this bias. In  \cite{GirolamiCalderhead11} a number of examples are provided illustrating the potential of such a scheme for challenging inference problems.

One can interpret the proposal mechanism of RMMALA as a local Gaussian approximation to the target density where the effective covariance matrix in RMMALA is the inverse of the metric tensor evaluated at the current position. Furthermore a simplified version of the RMMALA algorithm, termed sRMMALA, can also be derived by assuming a manifold with constant curvature thus cancelling the last term in Equation (\ref{eq:meanmMALA}) which depends on the Christoffel symbols. Whilst this is a step forward in that much information about the target density is now embedded in the proposal mechanism it is still driven by a random walk. The next approach to be taken goes beyond the direct and scaled random walk by defining proposals which follow the geodesic flows on the manifold of densities and thus presents a potentially really powerful scheme to explore posterior distributions.

\subsection{Riemann Manifold Hamiltonian Monte Carlo}
The Riemann manifold Hamiltonian Monte Carlo (RMHMC) method defines a Hamiltonian on the Riemann manifold of probability density functions by introducing the auxiliary variables $\pb\sim \mathcal{N}(\boldsymbol{0},\G(\ub))$ which are interpreted as the momentum at a particular position $\ub$ and by considering the negative log of the target density as a potential function. More formally the Hamiltonian defined on the Riemann manifold is
\begin{equation}
\Hc(\ub,\pb) = -\mathcal{L}(\ub) +\frac{1}{2}\log\left(2\pi|\G(\ub)|\right) + \frac{1}{2}\pb^T\G(\ub)^{-1}\pb
\end{equation}
where the terms $-\mathcal{L}(\ub) +\frac{1}{2}\log\left(2\pi|\G(\ub)|\right) $ and $\frac{1}{2}\pb^T\G(\ub)^{-1}\pb$ are the potential energy and kinetic energy terms respectively.
 and the dynamics given by Hamiltons equations are
\begin{eqnarray}
\frac{d \ub_k}{dt} &= \frac{\partial \Hc }{\partial \pb_k} = \left(\G(\ub)^{-1}\pb\right)_k \\
\frac{d\pb_k}{dt} &= -\frac{\partial \Hc }{\partial \ub_k} = \frac{\partial \mathcal{L}(\ub)}{\partial \ub_k} -\frac{1}{2}Tr\left[\G(\ub)^{-1}\frac{\partial \G(\ub)}{\partial \ub_k} \right] \nonumber\\
&+\frac{1}{2}\pb^T\G(\ub)^{-1}\frac{\partial \G(\ub)}{\partial \ub_k}\G(\ub)^{-1}\pb
\eqnlab{Hamilton}
\end{eqnarray}

These dynamics define geodesic flows at a particular energy level and as such make proposals which follow deterministically the most efficient path across the manifold from the current density to the proposed one. Simulating the Hamiltonian requires a time-reversible and volume preserving numerical integrator. For this purpose the Generalised Leapfrog algorithm can be employed and provides a deterministic proposal mechanism for simulating from the conditional distribution, i.e. $\ub^*|\pb \sim \pi(\ub^*|\pb)$. More details about the Generalised Leapfrog integrator can be found in  \cite{GirolamiCalderhead11}. To simulate a path (which turns out to be a local geodesic) across the manifold, the Leapfrog integrator is iterated $L$ times which along with the integration step size $\epsilon$ are parameters requiring tuning. Again due to the discrete integration errors on simulating the Hamiltonian in order to ensure convergence to the stationary distribution the Metropolis-Hastings acceptance ratio is applied.

The RMHMC method has been shown to be highly effective in sampling from posteriors induced by complex statistical models and offers the means to efficiently explore the hugely complex and high dimensional posteriors associated with PDE-constrained inverse problems.

\section{Low rank approximation of the Fisher information matrix}
\seclab{lowrank} As presented in Section \secref{HMC}, we use the
Fisher information matrix at the MAP point augmented with the Hessian
of the prior as the metric tensor in our HMC simulations. It is
therefore necessary to compute the augmented Fisher matrix and its
inverse. In~\cite{Bui-ThanhGhattas12a, Bui-ThanhGhattas12,
  Bui-ThanhGhattas12f}, we have shown that the Gauss-Newton Hessian of
the cost functional \eqnref{Cost}, also known as the data misfit, is a
compact operator, and that for smooth $\u$ its eigenvalues decay
exponentially to zero. Thus, the range space of the Gauss-Newton
Hessian is effectively finite-dimensional even before discretization,
i.e., it is independent of the mesh. In other words, the Fisher
information matrix admits accurate low rank approximations and the accuracy can be improved as desired by simply increasing the rank of the approximation. We shall
exploit this fact to compute the augmented Fisher information matrix
and its inverse efficiently. We start with the augmented Fisher
information matrix in $\R^N_\M$
\begin{align*}
\G := \M^{-1}\H + \alpha \Vb \boldsymbol{\Lambda}^{-2} \Vb^T\M = \alpha
\Vb\boldsymbol{\Lambda}^{-1}\LRp{\frac{1}{\alpha}\boldsymbol{\Lambda}\Vb^T\H\Vb\boldsymbol{\Lambda}
  + \I}\boldsymbol{\Lambda}^{-1}\Vb^{-1},
\end{align*}
where $\H$ is the Fisher information matrix obtained from
\eqnref{FisherF} by taking $\uo$ and $\ut$ as FEM basis functions.

Assume that $\H$ is compact (see, e.g., \cite{Bui-ThanhGhattas12a,
  Bui-ThanhGhattas12}), together with the fact that
$\boldsymbol{\Lambda}_{ii}$ decays to zero, we conclude that the
prior-preconditioned Fisher information matrix
\[
\Ht := \frac{1}{\alpha}\boldsymbol{\Lambda}\Vb^T\H\Vb\boldsymbol{\Lambda}
\] also has
eigenvalues decaying to zero. Therefore it is expected that the
eigenvalues of the prior-preconditioned matrix decays faster than
those of the original matrix $\H$. Indeed, the numerical results in
Section \secref{results} will confirm this observation.  It follows that (see,
e.g., \cite{FlathWilcoxAkcelikEtAl11, Bui-ThanhGhattasMartinEtAl13}
for similar decomposition) $\Ht$ admits a $r$-rank approximation of the
form
\[
\Ht = \frac{1}{\alpha}\boldsymbol{\Lambda}\Vb^T\H\Vb\boldsymbol{\Lambda} \approx
\Vb_r\S\Vb_r^T,
\]
where $\Vb_r$ and $\S$ (diagonal matrix) contain the first $r$
dominant eigenvectors and eigenvalues of $\Ht$, respectively. In this
work, similar to \cite{Bui-ThanhGhattasMartinEtAl13}, we use the
one-pass randomized algorithm in \cite{HalkoMartinssonTropp11} to
compute the low rank approximation.  Consequently, the augmented Fisher information matrix becomes
\[
\G \approx \alpha \Vb\boldsymbol{\Lambda}^{-1}\LRp{\Vb_r\S\Vb_r^T
  + \I}\boldsymbol{\Lambda}^{-1}\Vb^{-1},
\]
from which we obtain the inverse, by using the Woodbury formula \cite{GolubVan96},
\[
\G^{-1} \approx \frac{1}{\alpha} \Vb\boldsymbol{\Lambda}\LRp{\I - \Vb_r\D\Vb_r^T
  }\boldsymbol{\Lambda}\Vb^{-1},
\]
where $\D$ is a diagonal matrix with  $\D_{ii} = \S_{ii}/\LRp{\S_{ii}+1}$.

In the RMHMC method, we need to randomly draw the momentum variable as $\pb \sim
\GM{\mb{0}}{\G}$. If one considers
\[
\pb = \sqrt{\alpha} \Vb \Lambda^{-1}\bb + \sqrt{\alpha} \Vb\boldsymbol{\Lambda}^{-1}\Vb_r\S^{1/2}\cb,
\]
where $\bb_i,\cb_i \sim \GM{0}{1}$, then one can show , by inspection, that
$\pb$ is distributed by $\GM{\mb{0}}{\G}$.

\section{Numerical results}
\seclab{results}

For convenience, let us recall that the finite element (FEM) approximation of the
posterior is given as
\begin{equation}
\eqnlab{finalPosterior}
\post \sim \exp\LRp{-\half\snor{\db -
    \mc{G}\LRp{\u}}^2_\L} \times \exp\LRp{-\frac{\alpha}{2}\LRp{\ub-\ub_0}^T\M\Vb\boldsymbol{\Lambda}^2\Vb^T\M\LRp{\ub-\ub_0}},
\end{equation}
$\ub_0$ is the FEM nodal value of the prior mean function $\u_0$, $\M$
is the mass matrix, $\Vb$ the matrix of eigenvectors defined in
\eqnref{eigenmatrix}, $\boldsymbol{\Lambda}$ the diagonal matrix introduced in
\eqnref{eqnC}, $\L = \sigma^2\I$, $\db$ vector of observation data, and $\mc{G}\LRp{\u}$
the forward map given by the forward equation
\begin{align*}
-\Div\LRp{e^\u\Grad \w} &= 0 & \text{ in } \Omega \\
-e^\u\Grad \w \cdot \mb{n} &= Bi \,u &\text{ on } \partial \Omega \setminus \Gamma_{R}, \\
 -e^\u\Grad \w \cdot \mb{n} &= -1 &\text{ on } \Gamma_{R},
 \end{align*}
that is discretized by the $H^1$-conforming FEM method.

In this section, we study Riemann manifold Monte Carlo methods and
their variations to explore the posterior \eqnref{finalPosterior}. In
particular, we compare the performance of four methods: i) sRMMALA
obtained by ignoring the third derivative in RMMALA, ii) RMMALA, iii)
sRMHMC obtained by first computing the augmented Fisher metric tensor
at the MAP point and then using it as the constant metric tensor, iv)
RMHMC. For all methods, we form the augmented Fisher information
matrix exactly using \eqnref{FisherF} with $\uo$ and $\ut$ as finite
element basis vectors. For RMMALA and RMHMC we also need the
derivative of the metric tensor which is a third order tensor.  It can
be constructed exactly using \eqnref{tensor} with $\uo, \ut$ and
$\uth$ as finite element basis vectors. We also need extra work for
the RMHMC method since each Stormer-Verlet step requires an implicit
solve for both the first half of momentum and full position. For
inverse problems such as those considered in this paper, the fixed
point approach proposed in \cite{GirolamiCalderhead11} does not seem
to converge. We therefore have to resort to a full Newton
method. Since we explicitly construct the metric tensor and its
derivative, it is straightforward for us to develop the Newton
scheme. For all problems considered in this section, we have observed
that it takes at most five Newton iterations to converge.

Note that we limit ourselves in comparing these four methods in the
Riemannian manifold MCMC sampling family. Clearly, other methods are available,
we avoid ``unmatched comparison'' in terms of cost and the level of
exploiting the structure of the problem. Even in this limited family,
RMHMC is most expensive since it  requires not only third derivatives
but also implicit solves, but the ability in generating almost
independent samples is attractive and appealing as we shall show.

Though our proposed approach described in previous sections are valid
for any spatial dimension $d$, we restrict ourselves to a one
dimensional problem, i.e. $d = 1$, to clearly illustrate our points
and findings. In particular, we take $\Omega = \LRs{0,1}$, $\Gamma_R =
\LRc{1}$. We set $Bi = 0.1$ for all examples. As discussed in Section \secref{infiniteBayes},
for the Gaussian prior to be well-defined, we take $s = 0.6 > n/2 = 1/2$.

\subsection{Two-parameter examples}
\seclab{twoparameter}

We start our numerical experiments with two parameters. This will help
demonstrate various aspects of RMHMC which are otherwise too
computationally expensive for high dimensional problems. In
particular, two-parameter example allows us to compute the complete
third derivative tensor and perform the Newton method for each
Stormer-Verlet step. This in turn allows us to show the capability of the full
RMHMC over its simplified variants in tackling challenging posterior
densities in which the metric tensor changes rapidly.

In order to construct the case with two parameters we consider FEM with
one finite element. We assume that there is one observation point,
i.e. $K = 1$, and it is placed at the left boundary $x = 0$. In the
first example, we first take $s = 0.6$, $\sigma = 0.1$, and $\alpha =
0.1$. The posterior in this case is shown in Figure
\figref{posterior1}. We start by taking a time step of $0.02$ with $100$
Stormer-Verlet steps for both sRMHMC and RMHMC. The acceptance rate
for both methods is $1$. One would take a time step of $2$ for both
sRMMALA and RMMALA to be comparable with sRMHMC and RMHMC, but the
acceptance would be zero. Instead we take time step of $1$ so that the
acceptance rate is about $0.5$ for sRMMALA and $0.3$ for RMMALA. The
MAP point is chosen as the initial state for all the chains with
$5000$ sample excluding the first $100$ burn-ins. The result is shown
in Figure \figref{comparison1}.

\begin{figure}[h!tb]
  \subfigure[$s = 0.6$, $\alpha = 0.1$, and $\sigma = 0.1$]{
    \includegraphics[trim=1cm 6.5cm 2cm 7.3cm,clip=true,width=0.32\columnwidth]{ThermalFin1DNelements1al60alpha1m-1noise1m-1}

    \figlab{posterior1}
  }
  \subfigure[$s =
    0.6$, $\alpha = 1$, and $\sigma = 0.01$]{
    \includegraphics[trim=1cm 6.5cm 2cm 7.3cm,clip=true,width=0.32\columnwidth]{ThermalFin1DNelements1al60alpha1p0noise1m-2}

    \figlab{posterior2}
  }
  \subfigure[$s =
    0.6$, $\alpha = 0.1$, and $\sigma = 0.01$]{
    \includegraphics[trim=1cm 6.5cm 2cm 7.3cm,clip=true,width=0.32\columnwidth]{ThermalFin1DNelements1al60alpha1m-1noise1m-2}

    \figlab{posterior3}
}
  \caption{The contours of the posterior for three combinations of $s$
    (the prior smoothness), $\alpha$ (the ``amount'' of the prior), and $\sigma$ (the noise standard deviation).} \figlab{Posterior1}
\end{figure}

\begin{figure}[h!tb]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,width=0.97\columnwidth]{RMHMC_Thermal_ComparisonN5000Epsilon21al60alpha1m-1noise1m-1}
  \caption{Comparison of simRMMALA, RMMALA, simRMHMC, and RMHMC:
    chains with 5000 samples, burn-in of $100$, starting at the MAP
    point. In this example, $s = 0.6$, $\alpha = 0.1$, and $\sigma =
    0.1$. Time step is $\varepsilon = 1$ for simRMMALA and RMMALA, and
    $\varepsilon = 0.02$ with the number of time steps $L = 100$ for
    simRMHMC and RMHMC. In the left column: the exact synthetic
    solution is black, the sample mean is red, and the shaded region
    is the $95\%$ credibility region. In the middle column: blue is
    the trace plot for $\u_1$ while green is for $\u_2$. In the right
    column: red and black are the autocorrelation function for $\u_1$
    and $\u_2$, respectively.}  \figlab{comparison1}
\end{figure}

As can be seen, the RMHMC chain is the best in terms of mixing by
comparing the second column (the trace plot) and the third column (the
autocorrelation function ACF). Each RMHMC sample is almost uncorrelated to
the previous ones. The sRMHMC is the second best, but the samples are
strongly correlated compared to those of RMHMC, e.g. one uncorrelated
sample for every $40$. It is interesting to observe that the full RMMALA and sRMMALA have performance in terms of auto-correlation length that is qualitatively similar at least in the first $5000$
samples. This is due to the RMMALA schemes being driven by a single step random walk that cannot exploit fully the curvature information available to the geodesic flows of RMHMC, see rejoinder of \cite{GirolamiCalderhead11}.

 Note that it is not our goal to compare the
behavior of the chains when they converge. Rather we would like to
qualitatively study how fast the chains are well-mixed (mixing time). This is
important for large-scale problems governed by PDEs since
``unpredicted'' mixing time implies a lot of costly waste in PDE solves which
one must avoid. Though RMHMC is expensive in generating a
sample, the cost of generating an uncorrelated/independent sample
seems to be comparable to sRMHMC for this example. In fact, if we
measure the cost in terms of the number of PDE solves, the total number of PDE
solves for RMHMC is $42476480 $ while it is $ 1020002$ for sRMHMC, a
factor of $40$ more expensive. However, the cost in generating an
almost uncorrelated/independent sample is the same since sRMHMC generates one
uncorrelated sample out of $40$ while it is one out of one for RMHMC.

To see how each method distributes the samples we plot one for every
five samples in Figure \figref{comparison1Trajectory}. All methods
seem to explore the high probability density region very
well. This explains why the sample mean and the $95\%$ credibility
region are similar for all methods in the first column of Figure
\figref{comparison1}.

\begin{figure}[h!t!b!]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,width=0.97\columnwidth]{HelmholtzN5000BurnIn100Epsilon2Nsteps1001al60alpha1m-1noise1m-1Trajectory}
  \caption{Comparison MCMC trajectories (plot one for each five
    samples) among simRMMALA, RMMALA, simRMHMC, and RMHMC: chains with
    5000 samples, burn-in of $100$, starting at the MAP point. In this
    example, $s = 0.6$, $\alpha = 0.1$, and $\sigma = 0.1$. Time step
    is $\varepsilon = 1$ for simRMMALA and RMMALA, and $\varepsilon =
    0.02$ with the number of time steps $L = 100$ for simRMHMC and
    RMHMC.} \figlab{comparison1Trajectory}
\end{figure}

\begin{figure}[h!tb]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,width=0.97\columnwidth]{RMHMC_Thermal_ComparisonN5000Epsilon41al60alpha1p0noise1m-2}
    \caption{Comparison among simRMMALA, RMMALA, simRMHMC, and RMHMC:
    chains with 5000 samples, burn-in of $100$, starting at the MAP
    point. In this example, $s = 0.6$, $\alpha = 1$, and $\sigma =
    0.01$. Time step is $\varepsilon = 1$ for simRMMALA and RMMALA, and
    $\varepsilon = 0.04$ with the number of time steps $L = 100$ for
    simRMHMC and RMHMC. In the left column: the exact synthetic
    solution is black, the sample mean is red, and the shaded region
    is the $95\%$ credibility region. In the middle column: blue is
    the trace plot for $\u_1$ while green is for $\u_2$. In the right
    column: red and black are the autocorrelation function for $\u_1$
    and $\u_2$, respectively.}
  \figlab{comparison2}
\end{figure}

In the second example we consider the combination $s = 0.6$, $\sigma =
0.01$, and $\alpha = 1$ which leads to the posterior shown in Figure
\figref{posterior2}. For sRMHMC and RMHMC, we take time step
$\varepsilon= 0.04$ with $L = 100$ time steps, while it is $1$ for
both sRMMALA and RMMALA. Again, the acceptance rate is unity for both
sRMHMC and RMHMC while it is $0.65$ for sRMMALA and $0.55$ for RMMALA,
respectively.  The result for four methods is shown in Figure
\figref{comparison2}. As can be seen, this example seems to be easier
than the first one since even though the time step is larger, the
trace plot and the ACF looks better. It is interesting to observe that
sRMHMC is comparable with RMHMC (in fact the ACF seems to be a bit
better) for this example. As a result, RMHMC is more expensive than
sRMHMC for less challenging posterior in Figure
\figref{posterior2}. Here, by less challenging we mean that the
posterior is quite well approximated by a Gaussian at the MAP point,
e.g. the metric tensor is almost constant. This is true for the
posterior in Figure \figref{posterior2} in which the Gaussian prior
contribution is significant, i.e., $\alpha = 1$ instead of $\alpha =
0.1$. Conversely, the posterior is challenging if the metric tensor
changes rapidly. Similar to the first example, one also see that the
sample mean and the $95\%$ credibility region are almost the same for
all methods.

\begin{figure}[h!tb]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,width=0.97\columnwidth]{RMHMC_Thermal_ComparisonN5000Epsilon21al60alpha1m-1noise1m-2}
  \caption{Comparison among simRMMALA, RMMALA, simRMHMC, and RMHMC:
    chains with 5000 samples, burn-in of $100$, starting at the MAP
    point. In this example, $s = 0.6$, $\alpha = 0.1$, and $\sigma =
    0.01$. Time step is $\varepsilon = 0.7$ for simRMMALA and RMMALA, and
    $\varepsilon = 0.02$ with the number of time steps $L = 100$ for
    simRMHMC and RMHMC. In the left column: the exact synthetic
    solution is black, the sample mean is red, and the shaded region
    is the $95\%$ credibility region. In the middle column: blue is
    the trace plot for $\u_1$ while green is for $\u_2$. In the right
    column: red and black are the autocorrelation function for $\u_1$
    and $\u_2$, respectively.}
  \figlab{comparison3}
\end{figure}

In the third example we consider the combination $s = 0.6$, $\sigma =
0.01$, and $\alpha = 0.1$ which leads to a skinny posterior with a
long ridge as shown in Figure \figref{posterior3}. For sRMHMC and
RMHMC, we take time step $\varepsilon= 0.02$ with $L = 100$ time
steps, while it is $1$ for both sRMMALA and RMMALA. Again, the
acceptance rate is unity for both sRMHMC and RMHMC while it is $0.45$
for sRMMALA and RMMALA.  The result for four methods is shown in
Figure \figref{comparison3}. For this example, the RMHMC is more
desirable than sRMHMC since the cost to generate an
uncorrelated/independent sample is smaller for the former than the
latter. The reason is that the total number of PDEs solves for the
former is $40$ times more than the latter, but one out of
very sixty samples is uncorrelated/independent.

\subsection{Multi-parameter examples}
In this section we choose to discretize $\Omega = \LRs{0,1}$ with
$2^{10} = 1024$ elements, and hence the number of parameters is
$1025$. For all simulations in this section, we choose $s = 0.6$,
$\alpha = 10$, and $\sigma = 0.01$. For synthetic observations, we
take $K = 64$ observations at $x_j = (j-1)/2^6$, $j =
1,\hdots,K$. Clearly, using the full blown RMHMC is out of the
question since it is too expensive to construct the third derivative
tensor and Newton method for each Stomer-Verlet step. For that reason,
the sRMHMC becomes the viable choice.  As studied in Section
\secref{twoparameter}, though sRMHMC loses the ability to efficiently
sample from highly nonlinear posterior surfaces compared to the full
RMHMC it is much less expensive to generate a sample since it does not
require the derivative of the Fisher information matrix. In fact
sRMHMC requires to (approximately) compute the Fisher information at
the MAP point and then uses it as the fixed constant metric tensor
throughout all leap-frog steps for all samples. Clearly, the gradient
\eqnref{Gradient} has to be evaluated at each leap-frog step, but it
can be computed efficiently using the adjoint method presented in
Section \secref{adjoint}.

Nevertheless, constructing the exact Fisher information matrix
requires $2\times 1025$ PDEs solves. This is impractical if the
dimension of the finite element space increases, e.g. by refining the
mesh. Alternatively, due to the compactness of the Hessian of the
prior-preconditioned misfit as discussed in Section \secref{lowrank},
we can use the randomized singular value decomposition (RSVD)
technique \cite{HalkoMartinssonTropp11} to compute its low rank
approximations. Shown in Figure
\figref{eigenSpectrum} are the first 35 dominant eigenvalues of the
Fisher information matrix and its prior-preconditioned counterpart. We
also plot 20 approximate eigenvalues of the prior-preconditioned
Fisher information matrix obtained from the RSVD method. As can be
seen, the eigen spectrum of the prior-preconditioned Fisher
information matrix decays faster than that of the original one. This
is not surprising since the prior-preconditioned Fisher operator is a
composition of the prior covariance, a compact operator, and the
Fisher information operator, also a compact operator. The power of the
RSVD is clearly demonstrated as the RSVD result for the first 20
eigenvalues is very accurate.

\begin{figure}[h!tb]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,width=0.97\columnwidth]{GNspectrum1024al60alpha1000noise1}
  \caption{The eigen spectrum of the Fisher information matrix, the
    prior-preconditioned Fisher matrix, and the first 20 eigenvalues
    approximated using RSVD. Here, $s = 0.6$, $\alpha = 10$, and $\sigma
    = 0.01$.}  \figlab{eigenSpectrum}
\end{figure}

\begin{figure}[h!tb]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,height = 0.7\columnwidth, width=1\columnwidth]{FullLowRankFisher1024al60alpha1000noise1}
  \caption{MCMC results of three sRMHMC method with i) the low rank
    Gauss-Newton Hessian, ii) the exact Gauss-Newton Hessian, and iii)
    the full Hessian. In the figure are the empirical mean (red line),
    the exact distributed parameter used to generate the observation
    (black line), and $95\%$ credibility (shaded region).}
  \figlab{UQmultiparameter}
\end{figure}

Next, we perform the sRMHMC method using three different constant
metric tensors: i) the low rank Gauss-Newton Hessian, ii) the exact
Gauss-Newton Hessian, and iii) the full Hessian. For
each case, we start the Markov chain at the MAP point and compute
$5100$ samples, the first $100$ of which is then discarded as burn-ins. The
empirical mean (red line), the exact distributed parameter used to
generate the observation (black line), and $95\%$ credibility region
are shown in Figure \figref{UQmultiparameter}. As can be seen, the
results from the three methods are indistinguishable. The first sRMHMC
is the most appealing since it requires $2 \times 20 = 40$ PDE solves
to construct the low rank Fisher information while the others need
$2\times 1025$ PDE solves. For large-scale problems with computationally
expensive PDE solves, the first approach is the method of choice.

To further compare the three methods we record the trace plot of the
first two ($1$ and $2$) and the last two ($1024$ and $1025$) parameters in
Figure \figref{UQmultiparameterTrace}. As can be observed, the chains
from the three methods seem to be well-mixed and it is hard to see the
difference among them. We also plot the autocorrelation function for
these four parameters. Again, results for the three sRMHMC methods are
almost identical, namely, they generate almost uncorrelated samples. We
therefore conclude that low rank approach is the least computational
extensive, yet it maintains the attractive features of the original
RMHMC. As such, it is the most suitable method for large-scale
Bayesian inverse problems with costly PDE solves.

\begin{figure}[h!t!b!]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,height = 1.\columnwidth, width=1\columnwidth]{FullLowRankFisherTrace1024al60alpha1000noise1}
  \caption{MCMC results of three sRMHMC method with i) the low rank
    Gauss-Newton Hessian (left column), ii) the exact Gauss-Newton
    Hessian (middle column), and iii) the full Hessian at the MAP
    point (right column). In the figure are the trace plot of the
    first two ($1$ and $2$) and last two ($1024$ and $1025$) parameters}
  \figlab{UQmultiparameterTrace}
\end{figure}

\begin{figure}[h!t!b!]
    \includegraphics[trim=1cm 6.5cm 2cm 6.5cm,clip=true,height = 1.\columnwidth, width=1\columnwidth]{FullLowRankCorrelation1024al60alpha1000noise1}
  \caption{MCMC results of three sRMHMC method with i) the low rank
    Gauss-Newton Hessian (left column), ii) the exact Gauss-Newton
    Hessian (middle column), and iii) the full Hessian at the MAP
    point (right column). In the figure are the autocorrelation function plot of the
    first two ($1$ and $2$) and last two ($1024$ and $1025$) parameters}
  \figlab{UQmultiparameterCorrelation}
\end{figure}

\section{Conclusions and future work}
\seclab{conclusions} We have proposed the adoption of a computationally  inexpensive Riemann manifold
Hamiltonian Monte method to explore the posterior of large-scale
Bayesian inverse problems governed by PDEs in a highly efficient manner. We first adopt an infinite
dimensional Bayesian framework to guarantee that the inverse
formulation is well-defined. In particular, we postulate a Gaussian
prior measure on the parameter space and assume regularity for the
likelihood. This leads to a well-defined posterior distribution. Then,
we discretize the posterior using the standard finite element method
and a matrix transfer technique, and apply the RMHMC method on the
resulting discretized posterior in finite dimensional parameter
space. We present an adjoint technique to efficiently compute the
gradient, the Hessian, and the third derivative of the potential
function that are required in the RMHMC context.
This is at the expense of
solving a few extra PDEs: one for the gradient, two for a
Hessian-vector product, and four for the product of third order
derivative with a matrix.
For large-scale
problems, repeatedly computing the action of the Hessian and third order derivative
is too computationally expensive and this motivates us to design a
simplified RMHMC in which the Fisher information matrix is computed
once at the MAP point. We further reduce the effort by constructing
low rank approximation of the Fisher information using a randomized singular value decomposition technique. The effectiveness of
the proposed approach is demonstrated on a number of numerical results
up to $1025$ parameters in which the computational gain is about two
orders of magnitude while maintaining the quality of the original
RMHMC method in generating (almost) uncorrelated/independent samples.

For more challenging inverse problems with significant change of
metric tensor across the parameter space, we expect that sRMHMC with
constant metric tensor is inefficient. In that case, RMHMC seems to be a
better option, but it is too computational extensive for large-scale
problems. Ongoing work is to explore approximation of the RMHMC
methods in which we approximate the trace and the third derivative in
\eqnref{Hamilton} using adjoint and randomized techniques.

\section*{References}
\bibliographystyle{unsrt}

redpajama_set_name: RedPajamaArXiv


While I endeavour to ensure that the information on this website is correct and current, I do not warrant its completeness or accuracy; nor do I not commit to ensuring that the website remains available or that the material on this website is kept up-to-date.
To the maximum extent permitted by applicable law I exclude all representations, warranties and conditions (including, without limitation, the conditions implied by law of satisfactory quality, fitness for purpose and the use of reasonable care and skill).
My liability is limited and excluded to the maximum extent permitted under applicable law. I will not be liable for any direct, indirect or consequential loss or damage arising under these terms and conditions or in connection with my website, whether arising in tort, contract, or otherwise – including, without limitation, any loss of profit, contracts, business, goodwill, data, income, revenue or anticipated savings.
However, nothing in these terms and conditions shall exclude or limit my liability for fraud, for death or personal injury caused by my negligence, or for any other liability which cannot be excluded or limited under applicable law.
Access to certain areas of this website is restricted. I reserve the right to restrict access to some areas of this website, or indeed the entire website, at my discretion.
You must not use the website for any purposes related to marketing without my express written consent.
You must not use the website to copy, publish or send material which is illegal or unlawful, or material which could give rise to legal action under Canadian and other applicable laws. All material you copy, publish or send via the website must not be defamatory, obscene, indecent, hateful, discriminatory or inflammatory; such material must not infringe any person's intellectual property rights or rights of confidence, impinge upon any person's privacy, or constitute incitement to commit a crime; further, material must not be misleading, deceptive, sexually explicit, threatening, abusive, harassing or menacing.
In respect of all material that you may post on the website, you are granting worldwide, irrevocable, non-exclusive, royalty-free licence to use, reproduce, adapt, publish, translate and distribute such material in any media, together with the right to sub-licence such rights.
These terms and conditions may be amended from time-to-time. Please check this page regularly to ensure you are familiar with the current version.
These terms and conditions constitute the entire agreement between us, in relation to your use of my website.
This notice will be governed by and construed in accordance with Canadian law, and any disputes relating to this notice shall be subject to the exclusive jurisdiction of Canadian courts of law.
Due to unforeseen family circumstances, I am choosing to take a sudden and indefinite leave of absence from my counselling and coaching practices. I have very seriously considered the ramifications of this inescapable decision. Unfortunately, events out of my control have occurred and I need to support my loved ones during this time.
I apologize for the sudden nature of this news. I am going to provide names and contact information of other practitioners who I hold in high regard and other affordable therapeutic support alternatives (see list below).
At this point, I cannot predict my return. I am happy to arrange a time to speak with you by telephone (my cell: #778-859-1246) to respond to any questions you may have, referral options, brainstorm alternate practitioner possibilities and additional resources.
To assure you, my personal health is good. I do not want to cause worry. Simply put, out of respect for other people's privacy, this is a matter of unexpected life events occurring in a manner that requires my full presence and attention for my loved ones.
Thank you for being in my life – each and every one of you. You've taught me how to be a better human being and inspire every step I take.
*If you are experiencing a mental health emergency, please call 911.
For non-life-threatening mental health or addictions concerns, contact the Access and Assessment Centre (ACC) at VGH, available 24-7, 365 days a year. Phone, drop-in and outreach support available.
Counselling fees based on household income. Master's-level therapists. Program has a dedicated intake worker who can also refer to other counselling services or groups. Offices in Vancouver, Richmond, Burnaby, New Westminster and Surrey.
Professional counselling for LGBTQ individuals and couples.
Professional counselling for residents of the North Shore. Sliding Scale.
Reduced fee. Secular counselling services provided at the Vancouver Unitarian Centre by supervised volunteers with Master's degrees in psychology or psychology-related fields. Individual, couples and family counselling.
Adler Centre – Counselling Clinic – 604-742-1818 http://www.adlercentre.ca/clinic.html Sliding scale individual and couples counselling. Counselling provided by counselling psychology graduate students at the Adler Centre, supervised by an experienced clinician.
Free. Counselling provided by counselling psychology graduate students, supervised by a psychologist. Clinic runs from September to April.
Counselling services provided by doctoral student interns, supervised by registered psychologists. $10-$40 per hour.
Free counselling for the general public by counselling psychology graduate students, supervised by a psychologist.
Counselling for adults, children and youth provided by supervised graduate students in counselling psychology. Services at the Surrey clinic are free and at the Burnaby clinic are offered on a sliding scale.
Group therapy for depression, anxiety, stress. Individual mental health assessments. Referral through your family doctor, or via a walk-in clinic.
Individual, couple and family counselling available to Catholic and non-Catholic families. Sliding scale. Offices in Vancouver and Surrey.
Telephone coaching for people with mild-moderate depression with or without accompanying anxiety. Coaching available in English, Cantonese, French and Punjabi. Family doctor's referral required to access this program. Free.
Sliding scale counselling to the Jewish and non-Jewish community.
Counselling offered in Madarin, Cantonese, Korean and English with a focus on helping new immigrants of Chinese and other ethnic origins. Play and art therapy also available.
Individual, couple and family counselling using Bowen Family Systems Therapy. Lower- cost counselling provided by supervised interns.
By-donation ($20.00 min) feminist counselling services provided by supervised master's-level and intern narrative therapists.
Pay-by-donation counselling for individuals and families. Offices in Surrey and South Vancouver. Graduate-level counsellors and counselling interns.
Offers counselling and support services in the areas of youth and family, anger management, and sexual orientation/gender identity issues. Counselling provided by supervised interns completing their Master's of Counselling Degree. Free.
Surrey Youth Resources Centre, Community Counselling Clinic 604-592-6200 Free counselling provided by graduate student interns for individuals, youth and families (must be have a child or youth in the family to access services).
Individual therapy with certified counsellor ($85-$95 an hour) or intern ($50 hour). Group therapy also available.
Counselling for victims of crime, their immediate family members and witnesses. Contact for eligibility.
Free professional, multilingual counselling for any British Columbian struggling with their own or their family member's gambling.
Qmunity – Free Counselling Program – 604-684-5307 http://www.qmunity.ca/older-adults/resources/counselling/ Counselling for members of the GLBTQ communities.


Hazel Ledbetter, pictured left, with Ima Hogg at Round Top, circa 1970. Courtesy of The Dolph Briscoe Center for American History. Image available on the Internet and included in accordance with Title 17 U.S.C. Section 107.
Hazel's Lone Oak Cottage, named after Ledbetter, pictured after restoration at Winedale. Courtesy of The Dolph Briscoe Center for American History. Image available on the Internet and included in accordance with Title 17 U.S.C. Section 107.
Hazel Gawley Ledbetter's Grave. Image available on the Internet and included in accordance with Title 17 U.S.C. Section 107.
LEDBETTER, HAZEL GAWLEY (1902–1992). Hazel Gawley Ledbetter, daughter of Julian Hawkins and Irene (Harbuck) Gawley, was born in Houston, Texas, on September 25, 1902. She was educated in the public schools of Houston and at the College of Industrial Arts (now Texas Woman's University) in Denton. She married Dr. Paul Veal Ledbetter, grandson of Dr. A. A. Ledbetter, on June 10, 1924. They were the parents of three children.
She was active in civic affairs as an early member of River Oaks Garden Club and as president of the Women's Committee of the Houston Symphony Orchestra. During World War II she was in charge of the Camp and Hospital program at Camp Swift. As a longtime friend of Ima Hogg, she participated in the preparation of Bayou Bend and its collection to become a part of the Museum of Fine Arts, Houston.
Her interest in architecture and restoration led her to Round Top in 1959. She purchased several farms, houses, and old buildings and began a significant restoration program. In 1961 she purchased the property that is now known as the Winedale Historical Center. She offered the original stagecoach inn structure to the Museum of Fine Arts. When it could not be moved to Bayou Bend she sold the entire property to Ima Hogg who gave it to the University of Texas at Austin. In the late 1960s she sold all her remaining Round Top properties to Charles L. Bybee and his wife Faith.
She died in Houston on October 6, 1992, and is buried in Glenwood Cemetery. In 1993 she was awarded, posthumously, the Ima Hogg Award in recognition of her work at Round Top.
Fayette County Records. Houston Post, June 9, 10, 11, 1968. Hubert Roussel, The Houston Symphony Orchestra, 1913–1971 (Austin: University of Texas Press, 1972). Woman's Collection, Texas Women's University, Denton.
Handbook of Texas Online, Paul Gervais Bell, Jr., "LEDBETTER, HAZEL GAWLEY," accessed April 18, 2019, http://www.tshaonline.org/handbook/online/articles/flele.


text:
b. 25/03/1837 Glasgow, Scotland. d. 12/09/1875 Harberton, Devon
Francis Edward Henry Farquharson (1837-1875) was born on 25th March 1837 in Glasgow, Scotland. He was commissioned into the 42nd Highlanders in 1855 and served in the trenches at Sebastopol in the Crimean War. On 9th March 1858, the 20-year-old Lieutenant would be at the forefront of the attack on the Martiniere at Lucknow.
Attacking in two lines, the 42nd dashed down the grassy slope in order to minimise their losses. Following Colonel Alexander Cameron, Farquharson stormed a bastion mounting two guns, which he spiked and put out of commission. He then held this position until the following morning but suffered a severe wound.
He was invalided home and his VC citation was published on 21st June 1859. He received his medal on 4th January 1860 from the Queen Victoria at Buckingham Palace. Farquharson was promoted to Captain in 1862 and retired with the rank of major. He died at his home in Harberton, Devon on 12th September 1875, aged just 38. He was buried in St Peter's Churchyard, Harberton. His medals are held by the Black Watch Museum, Perth, Scotland.
LOCATION OF MEDAL: BLACK WATCH MUSEUM, PERTH, SCOTLAND.
BURIAL PLACE: ST PETER'S CHURCHYARD, HARBERTON, DEVON.
Francis Edward Henry Farquharson VC
Francis Farquharson's medals including his VC displayed at the Black Watch Museum, Perth, Scotland.
(Picture - Thomas Stewart).
Louis Desanges painting of Francis Farquharson VC


For over a year now, the wedding industry in Charlottesville has been anxiously awaiting a brand new venue to be built called "The Lodge". While new venues come all the time, this one was going to be built, owned, and operated by the same owners of Mount Ida which brought big shoes to fill (Mount Ida is already so stunning!). Shooting at Mount Ida last year, I saw some construction and wondered who would get to be the first videographer to grace the brand new venue for a wedding.
Enter Lauren and Michael. When Lauren and Michael inquired and within the second sentence were telling me about their boy cat named "Sassy", I new they would be a great couple to work with. Imagine my surprise when they told me they were getting married at the brand new Lodge! Lauren and Michael booked the venue before it was even done being built, and were the second couple ever to have a wedding there. Since the first wedding did not have a videographer, we were thrilled to know that we would be the first videographers to "break the ground".
The day couldn't had gone better, the day was beautifully orchestrated by Courtney and her staff from the Local, and I got to shoot alongside Meredith Sledge Photography, a photographer I have surprisingly never worked with, but hope to work with more in the future! Meredith's style and look is definitely something unique in the industry. I also always enjoy working with Andy Wilfong of Party Masters, because so few DJ's seem to be able to get me a clean signal out of their board, and not only is Andy able to do this but he is super helpful and goes out of his way to help us ensure a great audio capture.
Anyway, before I finish what is beginning to be a book: THE VOWS, GUYS! Lauren and Michael love each other so much, and while that is the case for many couples, they sure know how to put it into words and make every guest tear up while watching the ceremony (and me, later, in editing, possibly twice). This is one of my favorite wedding films because Lauren and Michaels love just engulfs every part of the day and is so easy and beautiful to capture. I ended going all the way back to the venue a week later (after I mastered my drone) and getting some shots of the venue. The Lodge is on thousands of acres far away from everything else, and I wanted the first shot of the film to reflect this amazing fact. I hope you enjoy this film! Grab a box of tissues though.