_Mhx<dZddlZddlmZdgZddZd Zed Zed Z d Z d Z dZ dZ ddZdZdZdZdZdS)zSparse block 1-norm estimator. N)aslinearoperator onenormestFct|}|jd|jdkrtd|jd}||krtjt|tj|}|j||fkr%tddt|jzt| d}|j|fkr%tddt|jztj |}t||} |dd|f} ||} nt||j||\} } } } } |s|r| f}|r|| fz }|r|| fz }|S| S)a Compute a lower bound of the 1-norm of a sparse array. Parameters ---------- A : ndarray or other linear operator A linear operator that can be transposed and that can produce matrix products. t : int, optional A positive parameter controlling the tradeoff between accuracy versus time and memory usage. Larger values take longer and use more memory but give more accurate output. itmax : int, optional Use at most this many iterations. compute_v : bool, optional Request a norm-maximizing linear operator input vector if True. compute_w : bool, optional Request a norm-maximizing linear operator output vector if True. Returns ------- est : float An underestimate of the 1-norm of the sparse array. v : ndarray, optional The vector such that ||Av||_1 == est*||v||_1. It can be thought of as an input to the linear operator that gives an output with particularly large norm. w : ndarray, optional The vector Av which has relatively large 1-norm. It can be thought of as an output of the linear operator that is relatively large in norm compared to the input. Notes ----- This is algorithm 2.4 of [1]. In [2] it is described as follows. "This algorithm typically requires the evaluation of about 4t matrix-vector products and almost invariably produces a norm estimate (which is, in fact, a lower bound on the norm) correct to within a factor 3." .. versionadded:: 0.13.0 References ---------- .. [1] Nicholas J. Higham and Francoise Tisseur (2000), "A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra." SIAM J. Matrix Anal. Appl. Vol. 21, No. 4, pp. 1185-1201. .. [2] Awad H. Al-Mohy and Nicholas J. Higham (2009), "A new scaling and squaring algorithm for the matrix exponential." SIAM J. Matrix Anal. Appl. Vol. 31, No. 3, pp. 970-989. Examples -------- >>> import numpy as np >>> from scipy.sparse import csc_array >>> from scipy.sparse.linalg import onenormest >>> A = csc_array([[1., 0., 0.], [5., 8., 2.], [0., -1., 0.]], dtype=float) >>> A.toarray() array([[ 1., 0., 0.], [ 5., 8., 2.], [ 0., -1., 0.]]) >>> onenormest(A) 9.0 >>> np.linalg.norm(A.toarray(), ord=1) 9.0 rz1expected the operator to act like a square matrixzinternal error: zunexpected shape axisN)rshape ValueErrornpasarraymatmatidentity Exceptionstrabssumargmaxelementary_vector_onenormest_coreH)Atitmax compute_v compute_wn A_explicit col_abs_sumsargmax_jvwestnmults nresamplesresults _/var/www/html/test/jupyter/venv/lib/python3.11/site-packages/scipy/sparse/linalg/_onenormest.pyrr sT AwqzQWQZLMMM  AAvvZ 0 3 3 : :2;q>> J JKK  1v % %.'#j.>*?*??AA A:***22  ! & &.'#l.@*A*AACC C9\** a * * qqq({ #8$(8ACE(J(J%Q6:I   qdNF   qdNF  cdfd}|S)z Decorator for an elementwise function, to apply it blockwise along first dimension, to avoid excessive memory usage in temporaries. c^|jdkr |S|d}tj|jdf|jddz|j}||d<~t |jdD] }|||z|||z<!|S)Nrrdtype)r r zerosr.range)xy0yj block_sizefuncs r(wrapperz%_blocked_elementwise..wrappers 71: " "477Na  n%%B!'!*!""5RXFFFAAkzkN:qwqz:>> < <$(D1Qz\>):$;$;!AjL.!!Hr))r6r7r5s` @r(_blocked_elementwiser9ys0 J       Nr)cn|}d||dk<|tj|z}|S)a9 This should do the right thing for both real and complex matrices. From Higham and Tisseur: "Everything in this section remains valid for complex matrices provided that sign(A) is redefined as the matrix (aij / |aij|) (and sign(0) = 1) transposes are replaced by conjugate transposes." rr)copyr r)XYs r( sign_round_upr>s4 AAa1fINA Hr)cRtjtj|dS)Nrr )r maxr)r<s r(_max_abs_axis1rAs 6"&))! $ $ $$r)c d}d}td|jd|D]?}tjtj||||zd}||}:||z }@|S)Nr+rr )r0r r rr)r<r5rr4r3s r(_sum_abs_axis0rDstJ A 1agaj* - - F26!Aa lN+,,1 5 5 5 9AA FAA Hr)cFtj|t}d||<|S)Nr-r)r r/float)rir"s r(rrs$ %   A AaD Hr)c|jdks|j|jkrtd|jd}tj|||kS)Nrz2expected conformant vectors with entries in {-1,1}r)ndimr r r dot)r"r#rs r(vectors_are_parallelrKsL v{{ag((MNNN  A 6!Q<<1 r)cb|jD]%tfd|jDsdS&dS)Nc38K|]}t|VdSNrK.0r#r"s r( z;every_col_of_X_is_parallel_to_a_col_of_Y..s.;;!'1--;;;;;;r)FT)Tany)r<r=r"s @r((every_col_of_X_is_parallel_to_a_col_of_YrUsL S;;;;qs;;;;; 55  4r)cj\}}dd|ftfdt|DrdS|"tfd|jDrdSdS)Nc3LK|]}tdd|fVdSrNrO)rQr4r<r"s r(rRz*column_needs_resampling..s: > > 1QQQT7 + + > > > > > >r)Tc38K|]}t|VdSrNrOrPs r(rRz*column_needs_resampling..s.77a#Aq))777777r)F)r rTr0rS)rGr<r=rrr"s ` @r(column_needs_resamplingrYs 7DAq !!!Q$A > > > > >U1XX > > >>>t} 777713777 7 7 4 5r)cztjdd|jddzdz |dd|f<dS)Nrrsizer)r randomrandintr )rGr<s r(resample_columnr_s>i1171:66q81rArdr@rJargsortrr)rATrA_linear_operatorAT_linear_operatorrr<g_prevh_prevkindr=gbest_jSZhr4s r(_algorithm_2_2rwsN)++)"--"A AA1uu9$$QAaC$99!;a?!!!QRR%qMAF F A ((C* J(//22 3 3 1  1  dddG !   J)0033 4 4 1   66!#a&&"&111f9qF|*L*LMM jmmDDbD!"1"% cFq 3 3A'3q622AaaadGG 66%fQi;; ? =>>>%fQi166 ? =>>> 661XX C C)!A$q ::C#$ABBBC QU*Z c6Mr)ct|}t|}|dkrtd|dkrtd|jd}||krtdd}d}tj||ft } |dkrjt d|D]} t| | t |D]7} t| | r%t| | |dz }t| | %8| t |z} tj dtj } d} tj ||ft } d}d} tj | | }|dz }t|}tj|}tj|}|| ks|dkr|dkr||}|dd|f}|dkr || kr| }n |} | }||krnt!|} ~t#| |rn|dkrIt |D]9} t| | |r&t| | |dz }t| | |&:~tj | | }|dz }t%|}~|dkrt|||krn.tj|ddd d|t)| z}~|dkrhtj|d|| rntj|| }tj||||f}t |D]}t3|||| dd|f< |d|tj|d|| }tj| |f} |dz }t3||}|||||fS) a Compute a lower bound of the 1-norm of a sparse array. Parameters ---------- A : ndarray or other linear operator A linear operator that can produce matrix products. AT : ndarray or other linear operator The transpose of A. t : int, optional A positive parameter controlling the tradeoff between accuracy versus time and memory usage. itmax : int, optional Use at most this many iterations. Returns ------- est : float An underestimate of the 1-norm of the sparse array. v : ndarray, optional The vector such that ||Av||_1 == est*||v||_1. It can be thought of as an input to the linear operator that gives an output with particularly large norm. w : ndarray, optional The vector Av which has relatively large 1-norm. It can be thought of as an output of the linear operator that is relatively large in norm compared to the input. nmults : int, optional The number of matrix products that were computed. nresamples : int, optional The number of times a parallel column was observed, necessitating a re-randomization of the column. Notes ----- This is algorithm 2.4. rz$at least two iterations are requiredrzat least one column is requiredrz't should be smaller than the order of Ar-NTrf)rr r r rhrFr0r_rYr/intprrrDr@rr>rUrArjlenr;isinall concatenater)rrkrrrlrmrr%r&r<rGind_histest_oldrtrprqr=magsr$rsind_bestr#S_oldrurvseenr4new_indr"s r(rrDsYR)++)"-- qyy?@@@1uu:;;;  AAvvBCCCFJ Ae$$$A 1uuq! " "A Aq ! ! ! !q  A)!Q// 1%%%a *!Q// qMAx)))HG !Qu%%%A A C= J(//22 3 3! a  fTll4 ==AFFAvvv;!!!V) A 66cWnnC  u99  !   3Au = =   q551XX $ $-aE::$#Aq)))!OJ.aE::$  J)0033 4 4!  1   66c!ff( ++  jmmDDbD!"21S]]?"2388:: q55ws2A2w))--// 73))D.#te*c$i!899Cq 3 3A'3q622AaaadGGbqb'273rr7H5556>8W"566 Q{=| !X&&A 1fj ((r))rrFFrN)__doc__numpyr scipy.sparse.linalgr__all__rr9r>rArDrrKrUrYr_rdrwrr8r)r(rs7000000 .kkkk\*     %%%           ===(((dddNO)O)O)O)O)r)