## Positive Matrix Factorization

Posted by Jeff Id on September 10, 2009

Well this was a good chunk of my weekend. It’s basically another method for solving the Antarctic reconstruction which is more of a brute force approach which has similarities to some optical calculations we do at my work. The method is conceptually simple and is performed in a few steps. The point of it was to attempt to avoid using PCA or other data reduction methods and come to a reasonable reconstruction.

It starts with something many have said before, AVHRR is surface skin temperature not surface air temperature and the methods of Steig et al fit ground information to the AVHRR data. Therefore they are reconstructing the surface skin temperature and not the surface air temperature. The result is presented as though they are the same and they well may be similar but I doubt the natural variance of temperature in an exposed rock hill or near meltpoint iceface is anywhere near as great as the variance in air temp. Anyway, I decided to first recalibrate the satellite data to ground data.

The steps are:

1 – Recalibrate AVHRR to better match ground data.

2 – Regress each gridcell of recalibrated AVHRR using the ground data employing iterative positive matrix factorization

My hope was to keep the calibration simple because complex versions simply confuse the situation. To do this the code collects from the area weighted offset reconstruction 64 complete infilled temperature records as created in this post Trends in the Antarctic

The linear difference between the slope of these 64 reconstructed temperatures and the satellite AVHRR temp was 0.13C/decade. Not terribly small. This linear difference was subtracted equally from all satellite series.

The code for calibration is here:

############################################################

############################################################

#recalibrate satellite data

#########################3#

#locate series from closest grid satellite data

satind=rep(NA,statnumber)for(i in 1:statnumber)

{

dist =circledist(surfcgrid[i,1],surfcgrid[i,2],sat.coord[,2],sat.coord[,1]) #dist is a vector of all stations

satind[i]=which(dist==min(dist))

}

satanom=calc.anom(all.avhrr)

satdata=ts(satanom[,satind],start=1982,deltat=1/12) #300×64

gnddata=ts(rec.area[,satind],start=1957,deltat=1/12) #600×64sd=ts(rowMeans(satdata),start=1982,deltat=1/12)

gd=ts(rowMeans(gnddata),start=1957,deltat=1/12)

gd=window(gd,start=1982)plot(sd,main=”Temperature Overlay Satellite and Ground \nAt 64 Ground Station Gridcells”)

lines(gd,col=”red”)

#lines(sd-gd,col=”green”,lwd=2)

abline(lsfit(time(sd),sd,NULL))

abline(lsfit(time(gd),gd,NULL),col=”red”)

sdcoef= round(coef(lsfit(time(sd),sd,NULL))[2]*10,3)

gdcoef= round(coef(lsfit(time(gd),gd,NULL))[2]*10,3)smartlegend(x=”left”,y=”bottom”,inset=0,fill=1:2,legend=c(paste(“Satellite=”,sdcoef,”C/Dec”),paste(“Ground”,gdcoef,”C/dec”)),col=c(“black”,”red”) ,cex=.7)

dslope=coef(lsfit(time(sd),sd-gd,NULL))[2]*10 #calculate deltaslope between ground and sat

correction=-(1:300)/120*dslope

correction=correction-mean(correction)

csatanom=satanom+correction #calculate recalibrated satellite data

csatdata=ts(csatanom[,satind],start=1982,deltat=1/12) #300×64

csd=ts(rowMeans(csatdata),start=1982,deltat=1/12)plot(csd,main=”Temperature Overlay Satellite and Ground \nAt 64 Ground Station Gridcells”)

lines(gd,col=”red”)

#lines(csd-gd,col=”green”,lwd=2)

abline(lsfit(time(csd),csd,NULL))

abline(lsfit(time(gd),gd,NULL),col=”red”)

sdcoef= round(coef(lsfit(time(csd),csd,NULL))[2]*10,3)

gdcoef= round(coef(lsfit(time(gd),gd,NULL))[2]*10,3)smartlegend(x=”left”,y=”bottom”,inset=0,fill=1:2,legend=c(paste(“Satellite=”,sdcoef,”C/Dec”),paste(“Ground”,gdcoef,”C/dec”)),col=c(“black”,”red”) ,cex=.7)

cor(csd,gd) #0.8294 recal sat data

cor(sd,gd) #0.8256 original sat data

The last two lines show that the corrected satellite data CSD has slightly better correlation after calibration. Since this step didn’t contain any offset corrections the two datasets after recalibration have the same slope and a noticeable offset.

Positive matrix factorization (PMF) was then used to find weighting coefficients to fit ground station data to the satellite data. The concept of positive matrix coefficient regression is likely very old but there was renewed interest in the 90’s which resulted in several papers on the topic. This method is my own creation derived from similar yet different methods used to calculate non-imaging optic lenses. The PMF moniker refers to the best fitting of a complex temperature profile in satellite data using ground thermometers with positive ONLY multiplication factors. The point is to get the best quality trend without the spurious correlation of inverted thermometers.

The algorithm is below. It contains several commented out lines which were used in debugging and experimentation but it runs as is.

############

# positive matrix factorization

############

pmf=function(gnd=gnd,sat=sat)

{

tgnd=window(gnd,start=1982)

cc=cor(tgnd,sat,use=”pairwise.complete.obs”)

#mask= cc >= 0gndm=tgnd#[,mask]

weights=rep(1,ncol(gndm))

delerror=1000000

curerror=1000000

cnt=100

steps=1

while (abs(delerror)>.0000001)

{

reco= rowMeans(t(t(gndm)*weights),na.rm=TRUE)

error=sat-reco

ee=cov(gndm,error,use=”pairwise.complete.obs”)

weights=weights+ ee^2*ee/abs(ee)*steps

weights[weights<0]=0

weights=as.vector(weights)/mean(weights)

olderror=curerror

curerror=(sum(error^2))

#print(curerror)

delerror=olderror-curerror

#print(c(delerror,cnt))

#print(c(delerror,cnt,steps))#steps=1/(abs(delerror)^0.5)*.05

if(cnt<5000)

{

steps=cnt/100

cnt=cnt+1

}else

{

print(c(delerror,cnt,steps))

if(abs(delerror)>.01)

{

print(steps)

return(rep(0,600))

}}

}

print(mean(weights))

print(curerror)

#fullrecon=rowMeans(t(t(gnd[,mask])*weights),na.rm=TRUE)fullrecon=rowMeans(t(t(gnd)*weights),na.rm=TRUE)

plot(reco,type=”l”)

lines(as.vector(sat),col=”red”)

abline(lsfit(1:300,reco,NULL))

abline(lsfit(1:300,sat,NULL),col=”red”)fullrecon

}

This method finds the error (difference between sat and ground) as a vector and calculates the covariance of all ground stations with the error quantity. This creates a weighted error vector ee with a value for each surface station. The value is positive for stations which have a bad correlation to the error and negative for stations which correlate well to the error. The ee vector values are squared keeping the sign intact and added to the weight vector. The weight vector is masked to insure that any values which have drifted negative are set to zero. The reconstruction of each individual satellite grid point is created by multiplying each surface station times a weight and finding the mean of all the stations. In most cases the mean weight averages to 1 simply by iteration, however sometimes the weights drifted up to 1.4 and down to 0.6 which since this is a temp scale has shifted away from thermometer calibrated degrees C. Therefore the weights were constrained further to always sum to 1 which turned out to make very little difference. Deltaerror is used to determine when convergence is reached and steps is designed to create more rapid convergence.

The algorithm takes about 4 hours to run through all the gridcells and I’ve run it partially and completely dozens of times now. The result of all this work is in Figure 3.

I was very disappointed with this result, especially after all the hours put in. It was supposed to look like the area weighted reconstruction shown in Figure 4.

Figure 4 is in my opinion the best representation of trend in the Antarctic over the last 50 years so I hope that new methods don’t end up distorting the trend. It’s clear that my method is blending high temp trends around the continent. There are several potential causes of this problem including spurious correlation and convergence to non-global minima. My own thought is that the noise level in the sat data is not allowing close enough matching to the ground stations and that results in weak regression results. Since the satellite data only exists post 1982 I looked at those comparisons as well.

Now we would expect a reasonable match from PMF trends to the satellite trends after running the algorithm. However, the AVHRR trends are far noisier than the ground station data and this effect has created strong trends simply through noise levels. Note the cooling section next to the extreme warming in the peninsula which doesn’t match the area weighted reconstruction. This is not a small problem. Another issue which I hadn’t noticed before becasue I thought Comiso had cleaned it up was the clear evidence of sea temperature contamination in shore pixels of the AVHRR. It’s obvious that the satellite is measuring some water temps in warm months sporadically creating pixels with strong different trends than their immediate neighbors. This is just more noise. Figure 6 is the 1982 onward PMF algorithm result.

This is visually a better match to the AVHRR data. I’m not terribly disturbed by the differences in trend distribution between the AVHRR and the ground data as the areas which are different are a better match to the area weighted reconstruciton.

My next step will likely be to do an improved and more serious calibration of the satellite data. Ryan did a good job of it removing offsets and trends in an older post, I may even go more severe than that correction. My hope is that the incorrect spreading of high trend stations across the continent will be more limited.

Finally, since the trends are clearly not quite right about the only thing we can get from this is the spatial distribution of the trends. The detail of the trends is not distorted in these images by PCA chiladni patterns so we can see which regions have warmed and cooled in Figure 3 and 6 although there is a positive bias to the trends in Fig 3. This can give an expected shape of distributions for the more sophisticated regularization alogrithms. Also, this is more evidence that the trend distribution patterns in Steig et al. were created by too few PC’s in the methods.

## Layman Lurker said

Regarding your figure #1, a key question for me has always been: Is this trend difference due to inherent differences in skin vs surface air temp, or is it due to bias introduced when the AVHRR data is processed in the cloud masking step.

I believe Antarctic clouds follow seasonal patterns. I wonder if there are seasonal differences in these trends.

## Jeff Id said

#1 Or is it created by poorly compensated transitions between satellite series?

My not know.

I think that if we subtract the difference or perhaps a heavily filtered difference out completely though it should make a reasonable correction that could improve fig 1.

## Jeff Id said

One of Ryan’s posts looked at the difference. It was pretty revealing in that several obvious steps were visible. Certainly we can reasonably remove at least those.

## Layman Lurker said

#2

Also, sat transition issues might have been amplified or further distorted by feeding false readings into cloud masking as well.

## Layman Lurker said

Here is Ryan’s difference plot from his “Antarctic Coup de Grace” post (which looks at ground stations only):

Note that your trend difference is even greater with AWS included in the calculation. It would be interesting to separate the AWS coastal locations from the inland locations to see if a) coastal AWS vs. ground stations differences with AVHRR are in the same ball park, and b) compare the effect of coastal and inland geography (where cloud effects are greatest) on the trend differences.

## Jeff Id said

#5 I verified this graph earlier this weekend using reconstructed gridcell data from the area weighted plot at each gridcell corresponding to the surface station.

I’m thinking now that a 6 to 24 month gaussian filter on the difference data subtracted from the entire AVHRR matrix might be a good calibration method.

## Ryan O said

Jeff,

Honestly, I think the 57-06 result looks pretty good. It seems like the PMF method may be a bit more influenced by grid cell variance than the fancier methods . . . it looks as if it coupled the coastal regions more than the weighted area, imputation, or RLS methods. But the general pattern is still the same.

Have you done verification stats for it?

## Jeff Id said

#7 No stats yet but they will be decent. I see each fit as they occur on a graph.

I’m running an improved version which looks like it will have almost the same result. It has much better fits though as the sat data is corrected by a filtered version of the sat/ground difference plot that Lurker referenced.

There are still a couple of hours left before it finishes though.

## Ryan O said

Cool.

BTW, with the correlation RLS, I’m getting CEs to entirely withheld stations of 0.68+. 0.58 was the max for covariance RLS. Writing a script to do a bunch of them so we can look at the graphs.

## Jeff Id said

Ryan,

I’ve made some more modifications. I’d like to see the quality of my results also. In Nic’s recent email he makes the point that the higher than expected trend may be due to something else. I’ve got to try subracting the filtered difference between satellite and ground stations at the 62 gridcells available right from the AVHRR. I see each gridcell fit as it happens during this reconstruction. It is a much much less black box calculation that way but it showed a visually huge improvement in quality of fit to each gridcell and limited some of the trend spreading.

Nic might be right on this portion of his point.

I’ve read every word of the back and forth discussion between both and don’t have a firm opinion yet. Scaling to SD of individual gridcells seems like a problem to me in that the noise level of the gridcell is probably a large chunk of the SD. I’ll have time this weekend as I’m at a trade show again.

## Nic L said

Jeff

Good work. I shall be interested to see the final results. I think that I follow your PMF method so far as the 64 grid cells containing ground stations go. One question. You say “The reconstruction of each individual satellite grid point is created by multiplying each surface station times a weight and finding the mean of all the stations.” I have obviously missed something, as it isn’t quite clear to me how you arrive at a full set of weights for all 5509 grid cells?

Regarding the email discussion between Ryan and myself, it may well be the case that scaling ground stations to the SD of the relevant individual grid cells doesn’t work well. It was just my attempt to arrive at something close to a RLS correlation solution, as I didn’t (and still don’t) see how an exact correlation solution can be achieved. I haven’t yet had a chance to study Ryan’s four recent emails, though.

## Geoff Sherrington said

Jeff,

You illustrate a perimeter problem where land meets sea. In some WIP at

http://landshape.org/enm/40-years-of-some-bom-australian-rural-temperature-data/

I’ve shown a perimeter problem for the last 40 years of Tav for the Australian mainland, for a too-small subset of rural stations. The coast sites have scarcely changed, inland sites have risen faster than an IPCC. The same effect appears on gridded animated GIFFs of temps around the globe, especially in SH around Africa, Aust and Sth America outside the tropics.

It gets more complex because Tmax and T mean have their own additional information content. Some yearly patterns converge, some diverge, some stay parallel at a given station. Thus the Tav is swung around by combining Tmax and Tmin.

These conclusions are from Australian Bureau data which have already been homogenised by them. If you are drawing data from secondary sources, there could likely be more adjustment that might not be real. They all derive from a primary BOM.

Australia has manned bases Casey, Mawson and Davis at the Antarctic perimeter since the IGY of the mid-1950s. These do not show much of a change over the last 40 years. For example, after a site shift at Davis, the 1971-2008 Tav from a linear fit shows a negative trend of 1.62 degrees a century, Casey from 1971-2008 gives positive 0.05 deg C a century and Mawson gives negative 0.21. Correlation coefficients between pairs of the 3 stations range from 0.5 to 0.9 on an annual basis. It is hard to reconcile this with a warming Antarctic, since the two outer bases are over 2,000 km apart.

Which brings me to Layman Lurker’s difference graph linked above. There are years like 2000 and 2005 where the difference is several degrees. One has to ask if this is noise or if it is real. I suspect that some is real, for an uncertain reason. That is the global hot peak of year 1998, which is not prominent or even apparent in the raw Antarctic data I have seen. The anomalous year of the surrounding decade at the three Aust bases was 1999, but it was decidedly colder than surrounding years. This admits the possibility that the Antarctic is not a uniform subset of the globe and so global assumptions should be applied cautiously. Specifically, the manned data might not be easily compatible with the satellite data. It then comes down to whether the differences are in precision (which, in an ideal world, would balance out) or bias, which would not.

## Ryan O said

Geoff, that graph is a

differencebetween satellite and ground temperatures at the ground station locations. So those 3 degree differences indicate a problem.But it’s not a graph of temperature anomalies – it’s a difference plot.

## Geoff Sherrington said

Ryan O, I understood it’s a difference plot.

The point that I’m trying to make is that some features of globaltemp records are present in ground readings (or 1.5m above) but not in satellite reconstructions. I used the 1988 global peak as an example. I can’t find it in the ground station data, though I have not looked at all sources. If a major event in a year can differ between Heaven and Earth, so it can over a century.