rlsqr

Share:

Description

R Version of LSQR routine written in Fortran by Paige and Saunders.

Usage

1
2
rlsqr(G=matrix(), u=vector(), wantse = 0, damp = 0, atol = 0,
 btol = 0, conlim = 0, itnlim = 100, nout = 0)

Arguments

G

Design Matrix

u

data vector

wantse

logical, if standard estimates are desired

damp

Damping parameter

atol

a Tolorance

btol

b Tolorance

conlim

con

itnlim

Iteration limit

nout

integer, output file (Not Available)

Details

This code is an R wrapper for performing the LSQR routine by Paige and Saunders. The LSQR program is a popular inversion program for solving the least squares problem, Ax=b.

Value

list:

x

solution vector

itn

number of iterations

xnorm

An estimate of the norm of the final solution vector x

rnorm

norm(rbar), function being minimized

arnorm

norm( Abar(transpose)*rbar), norm of the residual for the usual normal equations

acond

An estimate of cond(Abar)

anorm

An estimate of the Frobenius norm of Abar

istop

reason for stopping:0=exact solution, 1=okay, 2=damp is zero, 3=damp is nonzero, 4=condition problem, 5=iteration limit reached

Note

The fortran code has write statements that are not available in R-code according to the standard R documentation.

Author(s)

Jonathan M. Lees<jonathan.lees@unc.edu>, Keehoon Kim<keehoon@email.unc.edu>

References

Paige, C., and Saunders, M., LSQR: An algorithm for sparse linear equations and sparse least squares. Trans. Math. Software, 1982, 8, 43-71

See Also

sirt, art

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
## Not run: 


library(png)
library(PEIP)
shift=function(x,ns)
{
nsig=c(rep(0,ns),x[1:(length(x)-ns)])
return(nsig)
}

imagesc<-function(G, col=grey((1:99)/100), ... )
  {
    #########  plot an image after flipping and transposing
    ###   to match the matlab code
   d = dim(G)
   b = G[d[1]:1,]
   image(t(b), col=col, ...)
  }
img= readPNG('original_image.png')
im = matrix(as.vector(img) , 40000, 1)

load('rtest.RDATA')


#################################################
# Blurring the image (DATA == d)
#################################################

d0=G %*% im
d=d0


#################################################
# Solving damped least squares problem
# (    G     ) * x  =  ( d )
# ( damp * I )         ( 0 )
#################################################

# damp = lambda
lambda=0

# data vector (u)
u=d

# tolerance of iteration
atol=1.0e-6;btol=1.0e-6
#atol=.Machine$double.eps;btol=.Machine$double.eps

# taking condition number of G into consideration ( 0 means to ignore it )
conlim=0

# the number of iteration of LSQR
itnlim=1000;



#################################################
# Running LSQR
#################################################

# Executing LSQR
lx=rlsqr(G=G,u=u,damp=lambda,itnlim=itnlim,atol=atol,btol=btol,conlim=0)
# lx$x : model vector inverted from LSQR
# lx$xnorm : || x ||^2
# lx$rnorm : || Gx - d ||^2 (Norm of residual)
# lx$itn : the number of iteration to get the result


# Results of LSQR


layout(matrix(1:3, ncol=3))
imagesc(matrix(img,200,200),main='Original Image',axes=FALSE, xlab='',ylab='')
imagesc(matrix(d,200,200),main='Blurred Image',axes=FALSE,xlab='',ylab='')
imagesc(matrix(lx$x,200,200),
main=paste('LSQR solution Itr. = ',lx$itn),axes=FALSE,xlab='',ylab='')

## End(Not run)

Want to suggest features or report bugs for rdrr.io? Use the GitHub issue tracker.