automatically regularized nonnegative solutons for ill conditioned linear systems n.
Download
Skip this Video
Loading SlideShow in 5 Seconds..
AUTOMATICALLY REGULARIZED NONNEGATIVE SOLUTONS FOR ILL-CONDITIONED LINEAR SYSTEMS PowerPoint Presentation
Download Presentation
AUTOMATICALLY REGULARIZED NONNEGATIVE SOLUTONS FOR ILL-CONDITIONED LINEAR SYSTEMS

Loading in 2 Seconds...

play fullscreen
1 / 34

AUTOMATICALLY REGULARIZED NONNEGATIVE SOLUTONS FOR ILL-CONDITIONED LINEAR SYSTEMS - PowerPoint PPT Presentation


  • 228 Views
  • Uploaded on

AUTOMATICALLY REGULARIZED NONNEGATIVE SOLUTONS FOR ILL-CONDITIONED LINEAR SYSTEMS. RONDALL E. JONES Ktech Corporation, Albuquerque, NM 87185 rejones@sandia.gov Inverse Problems In Engineering Seminar 2006. Motivation. First look solution for sophisticated professionals

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
Download Presentation

PowerPoint Slideshow about 'AUTOMATICALLY REGULARIZED NONNEGATIVE SOLUTONS FOR ILL-CONDITIONED LINEAR SYSTEMS' - JasminFlorian


Download Now An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
automatically regularized nonnegative solutons for ill conditioned linear systems

AUTOMATICALLY REGULARIZED NONNEGATIVE SOLUTONS FOR ILL-CONDITIONED LINEAR SYSTEMS

RONDALL E. JONESKtech Corporation, Albuquerque, NM 87185rejones@sandia.gov

Inverse Problems In Engineering Seminar 2006

motivation
Motivation
  • First look solution for sophisticated professionals
  • Makes advanced methods available to non-expert
  • Enables adaptive regularization in automated contexts
historical background
Historical Background
  • Tiknohov regularization augments a linear system as follows:

Ax = b

λx = 0

  • Where the value to use for λ is problematic
  • The larger the value of λ the larger the residual and the smoother the solution
  • If the error in b is known approximately, the discrepancy principle says to adjust λuntil the RMS of the residual equals λ.
historical background1
Historical Background
  • Wahba (1970s)
  • Generalized Cross Validation
  • Heuristic
  • Not reliable enough
historical background2
Historical Background
  • Jones (1980s)
  • Picard condition plus a simple error estimate and use of discrepancy principle
  • Fundamental analysis
  • Not sufficiently developed for general use
historical background3
Historical Background
  • Hansen (1990s)
  • Automated L-curve analysis
  • Heuristic
  • Quite good; but software not widely available
historical background4
Historical Background
  • O-Leary (2000s)
  • Based on certain assumptions, minimized the “distance” of the regularized solution to the true solution
  • Fundamental analysis; unbiased, which turns out not necessarily an advantage
historical background5
Historical Background
  • Jones (2005+)
  • More sophisticated Picard analysis; more robust error analysis
  • Fundamental analysis
  • Software package with realistic constraints for general usage
the picard condition
The Picard Condition
  • Ax = b (with A any shape)
  • A=USVT (using “compact version of SVD”)
  • USVT x = b
  • VT x = S+UTb == t, a vector
  • Si 0 in ill-conditioned case; (UTb)i epsilon
  • x = Vt
  • The Picard condition says that the magnitudes of the ti should generally decrease
  • If they instead grow after some point, then noise dominates the bi and ti and the solution will be bad
typical t vector for highly ill conditioned problem hansen s ilaplace
Typical t vector for Highly Ill-conditioned Problem(Hansen’s ilaplace)
usable rank determination
Usable Rank Determination
  • So how to determine when to cut off more elements of t as noise dominated?
  • Typical approach is to compute a moving average of the (absolute values of the) ti centered at each ti
  • But the average will then have an odd number of terms and will echo any oscillation
usable rank determination solution
Usable Rank Determination: Solution
  • Use an even number of terms in the averages
  • Find the minimum of that series, avgm
  • See if the average ever gets “significantly larger”
  • If not, consider it stable.
  • If so, find the first average, avgk, which is 3x larger than the minimum
  • Look inside that avgk’s span to find the largest ti and set ur=i-1
  • Backtrack ur as far as the start of the minimum average as long as the ti are backing down
  • This is obviously heuristic
our algorithm
Our Algorithm
  • Compute the SVD and the usable rank, ur, as described
  • Compute an estimate of σ, the error in the right hand side, as follows
  • Use the discrepancy principle to pick λ
  • Apply constraints
  • Only the determination of ur is heuristic
simple error estimation
Simple Error Estimation
  • Ax = b (with A any shape)
  • USVT x = b
  • VT x = S+UTb == t
  • {VT x = t}1ur usable portion
  • {VT x = t}ur+1m noise-dominated portion
  • RMS(t(ur+1:m)) is a decent estimate of the error (the standard deviation of the elements of b)
simple error estimation1
Simple Error Estimation
  • Problem: there may be as few as ONE element in RMS(t(ur+1:m))
  • Too crude… and there is a better way
better error estimation
Better Error Estimation
  • Ax = b original
  • {VT x = t}1ur orthogonalized usable part
  • Subtract from each row of (1) its projection on every row of (2) leaving Ax = b where b has essentially no image of the true b… just image of the noise
  • RMS(b) is almost an estimate of error in b…
  • But the bi are biased estimates because VT was derived from A. So we must adjust for lost degrees of freedom.
better error estimation1
Better Error Estimation
  • To correct for bias, use typical adjustment for lost degrees of freedom

i=m

  • σ2 = (1/(m-ur)) sum bi2

i=1

  • The identical formula can be obtained by replacing the orthogonalized usable part with a row-wise orthogonalized decomposition of A (with b carried along)
  • The correctness of this formula is rapidly seen in practice
our algorithm1
Our Algorithm
  • Compute the SVD and the usable rank, ur, as described
  • Compute the estimate of σ as just described
  • Use the discrepancy principle to pick λ
  • Apply constraints
non negativity
Non-Negativity
  • Traditional approaches exist… but rarely for the ill-conditioned case
  • We researched several approaches
  • We need some degree of simplicity because non-negativity is just the first phase of constraint we intend
non negativity1
Non-Negativity
  • I experimented with non-traditional approaches. For example…

a11 a12 a13 a14 a15 x1 b1

a21 a22 a23 a24 a25 x2 b2

a31 a32 a33 a34 a35 x x3 = b3

a41 a42 a43 a44 a45 x4 b4

a51 a52 a53 a54 a55 x5 b5

  • x4 can be made zero by subtracting from b its projection onto column 4 of A.
  • To see this, note X = (ATA)-1ATb. The 4th component of ATb will be zero after the projection
non negativity2
Non-Negativity
  • The problem with such techniques is that there is no control over the increase in problem residual. And the residual typically increases dramatically more than with column zeroing.
non negativity3
Non-Negativity
  • Simple setting of a variable (and its column) to zero was compelling
  • But in what sequence?
    • Experimentation suggests that the most effective variable to set to zero is the one which is hardest to adjust to zero by adjusting the RHS… which is usually the most negative value
  • We recognize the existence of much more sophisticated non-negativity algorithms. But we question the added value in our ill-conditioned context with other constraints to be added.
approach
Approach
  • Perform autoreg algorithm
  • Iteratively,

- Zero most negative xi (i.e., its column of A)

    • (in massive problems, zero several)
    • Re-solve regularized problem using the λ determined originally
other constraints equality
Other Constraints - Equality
  • Equality constraints can be implemented by subtracting from each regular equation its projection onto each of the equality constraint equations.
  • Then the reduced problem can be solved as before.
  • The solutions to parts 1 and 2 are added.
  • Some theoretical questions about this process do arise.
  • What if the equality constraints are “bad”? That is, redundant, conflicting, or ill-conditioned?
  • Our software handles all these cases comparatively easily.
other constraints inequality available in rmatrix
Other Constraints – Inequality(available in RMatrix)
  • Inequality constraints can be implemented by iteratively selecting certain of the inequality constraints to be advanced to equality status.
  • Reasoning as for non-negativity, the natural choice of next inequality constraint to advance seems to be the most-violated constraint.
  • The algorithm is to do an initial solution to pick lambda, then iteratively pick an inequality constraint to advance to equality status, and resolve using the equality constraint algorithm.
  • Some theoretical questions about this process arise. These have not been adequately studied.
  • What if the inequality constraints are “bad”? That is, redundant, conflicting, or ill-conditioned?
  • Our software handles all these cases. Basically, the first constraint selected overrides others which conflict or are redundant.
beyond tikhonov
Beyond Tikhonov
  • Tikhonov solves

Ax = b

λIx = 0

  • But sometimes the engineer wants a matrix other than I, such as a matrix representing the solution’s first, second, or third derivative rather than the solution itself.

Ax = b

λLx = 0 where L1(1)= (-1 1 0 0 0 0 …)

or L2(1) =( 1 -2 1 0 0 0 …)

or L3(1) =(-1 3 -3 1 0 0 … )

beyond tikhonov1
Beyond Tikhonov
  • The efficiency of the SVD when searching for the correct lambda (only one decomposition for all lambdas tried for one solution) is not available for non-identity L matrices.
  • So the computation is relatively expensive.
  • Non-negativity can be achieved for these cases in the same manner as Tikhonov. But now the computation is quite expensive because many SVDs (or QRs or whatever) must be done for each newly zeroed column.
  • The offered software provides all these features.
  • Second-derivative L matrix looks very useful.
caveats
Caveats
  • No automatic method can ever be perfect
    • The “bad” oscillations could be right!
  • Our heuristics have been tuned on a variety of problems
  • Always check the answer for suitability
software
Software
  • Autoreg3 is in beta test
  • Packaged as a single *.h file
  • Email the author for a copy
  • Matlab versions are available
  • RMatrix version is in commercial evaluation
  • These are NOT public domain
  • Wed site will be coming