Critérios de parada para solucionadores lineares iterativos aplicados a sistemas quase singulares

Considere com quase singular, o que significa que há um autovalor de que é muito pequeno. O critério de parada usual de um método iterativo é baseado no residual e diz respeito às iterações podem parar quando com o número da iteração. Mas no caso que estamos considerando, pode haver um grande erro vivendo no espaço próprio associado ao pequeno autovalor que fornece um pequeno valor residual . Suponha que o residual inicial seja grande, então pode acontecer que paremos em $Ax=b$ $A$ $\lambda_0$ $A$ $r_n:=b-Ax_n$ $\|r_n\|/\|r_0\|<tol$ $n$ $v$ $\lambda_0$ $Av=\lambda_0v$ $r_0$ $\|r_n\|/\|r_0\|<tol$ mas o erro $x_n-x$ ainda é grande. Qual é o melhor indicador de erro neste caso? É $\|x_{n}-x_{n-1}\|$ um bom candidato?

linear-algebra Hui Zhang
fonte

Você pode pensar na sua definição de "quase singular". A matriz

I \cdot ϵ

$I \cdot \epsilon$ (com

ϵ ≪ 1

$\epsilon\ll 1$ e

I

$I$ a matriz de identidade) tem um valor próprio muito pequeno, mas está tão longe de ser singular quanto qualquer matriz poderia estar.

David Ketcheson

Além disso,

| | r_{n} / r_{0} | |

$||r_n/r_0||$ parece a notação errada.

| | r_{n} | | / | | r_{0} | |

$||r_n||/||r_0||$ é mais típico, não?

Bill Barth

Sim, você está certo, Bill! Eu vou corrigir esse erro.

Hui Zhang

E quanto a? e qual é o seu algoritmo, exatamente?

‖ b - A x ‖ / ‖ b ‖

$\| b - Ax \| / \| b \|$

shuhalo 6/01/12

Adendo: Eu acho que o artigo a seguir aborda praticamente os sistemas mal condicionados com os quais você se preocupa, pelo menos se você estiver usando CG: Axelson, Kaporin: Estimativa da norma de erro e critérios de parada em iterações de gradiente conjugado pré-condicionadas. Você pode usar o seguinte

código

Respostas:

Por favor, nunca use a diferença entre iterações sucessivas para definir um critério de parada. Isso diagnostica estagnação por convergência. A maioria das iterações não simétricas da matriz não é monótona e mesmo o GMRES na aritmética exata sem reinicializações pode estagnar por um número arbitrário de iterações (até a dimensão da matriz) antes de convergir repentinamente. Veja exemplos em Nachtigal, Reddy e Trefethen (1993) .

Uma maneira melhor de definir convergência

Geralmente, estamos interessados na precisão da nossa solução mais do que no tamanho do resíduo. Especificamente, gostaríamos de garantir que a diferença entre uma solução aproximada e a solução exata satisfaça para alguns especificados pelo usuário . Acontece que é possível conseguir isso encontrando um tal que que é o menor valor singular de , devido a $x_n$ $x$

| x_{n} - x | < c

$|x_n - x| < c$

c

$c$

x_{n}

$x_n$

| A x_{n} - b | < c ϵ

$|A x_n - b| < c\epsilon$

ϵ

$\epsilon$

A

$A$

\begin{aligned} | x_{n} - x | & = | A^{- 1} A (x_{n} - x) | \\ \leq \frac{1}{ϵ} | A x_{n} - A x | \\ = \frac{1}{ϵ} | A x_{n} - b | \\ < \frac{1}{ϵ} \cdot c ϵ = c \end{aligned}

$\begin{align} |x_n - x| &= |A^{-1} A (x_n - x)| \\ & \le \frac 1 \epsilon |A x_n - A x| \\ & = \frac 1 \epsilon |A x_n - b| \\ & < \frac 1 \epsilon \cdot c \epsilon = c \end{align}$

onde usamos que é o maior valor singular de (segunda linha) e que resolve exatamente (terceira linha). $1/\epsilon$ $A^{-1}$ $x$ $A x = b$

Estimando o menor valor singular $\epsilon$

Uma estimativa precisa do menor valor singular geralmente não está disponível diretamente no problema, mas pode ser estimada como um subproduto de um gradiente conjugado ou iteração GMRES. Observe que, embora as estimativas dos maiores valores próprios e valores singulares sejam geralmente muito boas depois de apenas algumas iterações, uma estimativa precisa do menor valor próprio / próprio / singular geralmente é obtida apenas quando a convergência é alcançada. Antes da convergência, a estimativa geralmente será significativamente maior que o valor real. Isso sugere que você deve realmente resolver as equações antes de definir a tolerância correta . Uma tolerância de convergência automática que requer uma precisão fornecida pelo usuário $\epsilon$ $c\epsilon$ $c$ para a solução e estimativas, o menor valor singular com o estado atual do método de Krylov pode convergir muito cedo, porque a estimativa de era muito maior que o valor real. $\epsilon$ $\epsilon$

Notas

A discussão acima também funciona com substituído pelo operador pré-condicionado à esquerda e pelo residual pré-condicionado ou com o operador pré-condicionado à direita e o erro . Se for um bom pré-condicionador, o operador pré-condicionado estará bem condicionado. Para o pré-condicionamento esquerdo, isso significa que o residual pré-condicionado pode ser reduzido, mas o residual real pode não ser. Para o pré-condicionamento correto,é facilmente reduzido, mas o verdadeiro erro $A$ $P^{-1}A$ $P^{-1} (A x^n - b)$ $A P^{-1}$ $P (x_n - x)$ $P^{-1}$ $|P(x_n - x)|$ $|x_n-x|$ talvez não seja. Isso explica por que o pré-condicionamento à esquerda é melhor para reduzir o erro, enquanto o pré-condicionamento à direita é melhor para reduzir o residual (e para depurar pré-condicionadores instáveis).
Veja esta resposta para mais discussões sobre normas minimizadas pelo GMRES e CG.
As estimativas de valores extremos singulares podem ser monitoradas usando-se -ksp_monitor_singular_valuecom qualquer programa PETSc. Consulte KSPComputeExtremeSingularValues () para calcular valores singulares do código.
Ao usar o GMRES para estimar valores singulares, é crucial que as reinicializações não sejam usadas (por exemplo, -ksp_gmres_restart 1000no PETSc).

Jed Brown
fonte

'' também funciona com A substituído por um operador pré-condicionado '' - No entanto, ele se aplica somente ao

residual pré-condicionado se

for usado, resp. ao erro pré-condicionado

for usado.

P^{- 1} r

$P^{-1}r$

P^{- 1} A

$P^{-1}A$

P^{- 1} δ x

$P^{-1}\delta x$

A P^{- 1}

$AP^{-1}$

Arnold Neumaier

Bom ponto, eu editei minha resposta. Observe que a caixa pré-condicionada à direita permite controlar , desenrolando o pré-condicionador (aplicando

) normalmente amplifica os modos de baixa energia no erro.

P δ x

$P\delta x$

P^{- 1}

$P^{-1}$

precisa

$Ax=b$ $\min ||Ax-b||_2$ onde está muito mal condicionado (ou seja, a razão entre o primeiro e o último valor singular é grande). $A$ $\sigma_1/\sigma_n$

Aqui, temos vários métodos para escolher o critério de parada e, para um método iterativo, eu recomendaria o critério da curva L, uma vez que envolve apenas quantidades que já estão disponíveis (AVISO LEGAL: Meu orientador foi pioneiro neste método, por isso estou definitivamente inclinado a isto). Eu usei isso com sucesso em um método iterativo.

A idéia é monitorar a norma residual e a norma de solução , onde é o 'é iterado. À medida que você itera, isso começa a desenhar a forma de um L em um gráfico de loglog (rho, eta), e o ponto no canto desse L é a escolha ideal. $\rho_k=||Ax_k-b||_2$ $\eta_k=||x_k||_2$ $x_k$ $k$

Isso permite que você implemente um critério no qual você fica de olho quando passa pela esquina (ou seja, observa o gradiente de ) e, em seguida, escolhe a iteração que estava localizada na esquina. $(\rho_k,\eta_k)$

A maneira como eu fiz isso envolveu o armazenamento das últimas 20 iterações, e se o gradiente foi maior que um limite para 20 iterações sucessivas, eu sabia que estava na parte vertical da curva e que havia passado a esquina. Então, tomei a primeira iteração da minha matriz (ou seja, a de 20 iterações atrás) como minha solução. $abs(\frac{\log(\eta_k)-\log(\eta_{k-1})}{\log(\rho_k)-\log(\rho_{k-1})})$

Também existem métodos mais detalhados para encontrar a esquina, e eles funcionam melhor, mas exigem o armazenamento de um número significativo de iterações. Brinque um pouco com isso. Se você estiver no matlab, poderá usar a caixa de ferramentas Regularization Tools, que implementa parte disso (especificamente a função "corner" é aplicável).

Observe que essa abordagem é particularmente adequada para problemas de grande escala, pois o tempo de computação extra envolvido é minúsculo.

OscarB
fonte

Muito obrigado! Então, no loglog (rho, eta), começamos a partir da direita da curva L e terminamos no topo de L, não é? Eu simplesmente não conheço o princípio por trás desse critério. Você pode explicar por que sempre se comporta como uma curva em L e por que escolhemos o canto?

Hui Zhang

Vocé é bem vindo. Para um método iterativo, começamos da direita e terminamos sempre no topo. Ele se comporta como um L devido ao ruído no problema - a parte vertical acontece em

, onde

é o vetor de ruído

| | A x - b | |_{2} = | | e | |_{2}

$||Ax-b||_2=||e||_2$

e

$e$

b_{e x a c t} = b + e

$b_{exact}=b+e$ . Para mais análises, ver Hansen, PC, & O'Leary, DP (1993). O uso da curva L na regularização de problemas discretos e mal colocados. Jornal SIAM de Computação Científica, 14. Observe que eu apenas fiz uma pequena atualização no post.

OscarB

@HuiZhang: nem sempre é um L. Se a regularização for ambígua, pode ser um L duplo, levando a dois candidatos à solução, um com o desempenho bruto melhor resolvido e o outro com alguns detalhes melhor resolvidos. (E, claro, mor formas ecomplex pode aparecer.)

Arnold Neumaier

A curva L se aplica a problemas mal condicionados, onde deve haver uma solução única? Ou seja, estou interessado nos problemas Ax = b, onde b é conhecido "exatamente" e A é quase singular, mas ainda tecnicamente invertível. Parece-me que, se você usar algo como GMRES, a norma de seu palpite atual de x não mudará muito ao longo do tempo, especialmente após a primeira, no entanto, muitas iterações. Parece-me que a parte vertical da curva L ocorre porque não há solução única / válida em um problema mal colocado; esta característica vertical estaria presente em todos os problemas mal condicionados?

Nkeguy 06/06/19

Em um ponto, você alcançará uma linha vertical, normalmente porque os erros numéricos no método de solução resultam em || Ax-b || não diminuindo. No entanto, você está certo de que, em problemas sem ruído, a curva nem sempre se parece com um L, o que significa que você normalmente tem alguns cantos para escolher e escolher um sobre o outro pode ser difícil. Acredito que o artigo que mencionei no meu comentário acima discuta brevemente os cenários sem ruído.

OscarB 7/06/16

Critérios de parada para solucionadores lineares iterativos aplicados a sistemas quase singulares

Respostas:

Uma maneira melhor de definir convergência

Estimando o menor valor singularϵϵ\epsilon

Notas

Estimando o menor valor singular $\epsilon$