copyA = np.random.rand(5, 3)   # (m,n)

# SVD completa
U, s, Vh = np.linalg.svd(A, full_matrices=True)
# U:  (m,m) ortogonal
# s:  (k,) valores singulares k=min(m,n), orden desc
# Vh: (n,n) conjugado transpuesto de V

# SVD reducida (economica) — mas eficiente
U, s, Vh = np.linalg.svd(A, full_matrices=False)
# U:(m,k), s:(k,), Vh:(k,n)

# Reconstruir A exactamente
A_rec = U @ np.diag(s) @ Vh

# Aproximacion de rango r (compresion)
r = 2
A_r = U[:, :r] @ np.diag(s[:r]) @ Vh[:r, :]

# Solo valores singulares (mas rapido)
s_only = np.linalg.svdvals(A)

copyA = np.random.rand(5, 3)

# QR reducida (default) — la mas usada
Q, R = np.linalg.qr(A, mode='reduced')
# Q: (m,k), R: (k,n) donde k=min(m,n)

# QR completa
Q, R = np.linalg.qr(A, mode='complete')
# Q: (m,m) ortogonal, R: (m,n) triang. sup.

# Verificar ortogonalidad
np.allclose(Q.T @ Q, np.eye(Q.shape[1]))  # True

# Resolver Ax=b con QR (mas estable que inv)
b = np.random.rand(5)
Q, R = np.linalg.qr(A)
x = np.linalg.solve(R, Q.T @ b)

copy# Crear matriz SPD garantizada
M = np.random.rand(4, 4)
A = M.T @ M + np.eye(4)   # siempre SPD

# Cholesky lower (default)
L = np.linalg.cholesky(A)        # A = L @ L.T

# Cholesky upper triangular
U = np.linalg.cholesky(A, upper=True)  # A = U.T @ U

# Verificar
np.allclose(L @ L.T, A)  # True

# Resolver Ax=b (2x mas rapido que LU para SPD)
b = np.random.rand(4)
y = np.linalg.solve(L, b)      # forward substitution
x = np.linalg.solve(L.T, y)   # back substitution

copyA = np.array([[3,1],[0,2]])

vals, vecs = np.linalg.eig(A)
# vals: (n,) eigenvalores (pueden ser complejos)
# vecs: (n,n) columnas son eigenvectores normalizados
# vecs[:,i] es el eigenvector de vals[i]

# Verificar Av = lambda*v
i = 0
np.allclose(A @ vecs[:,i], vals[i] * vecs[:,i])  # True

# Diagonalizacion: A = V @ diag(lambda) @ V_inv
V = vecs
V_inv = np.linalg.inv(V)
np.allclose(V @ np.diag(vals) @ V_inv, A)  # True

# Ordenar por eigenvalor descendente
idx = np.argsort(np.abs(vals))[::-1]
vals_sorted = vals[idx]
vecs_sorted = vecs[:, idx]

copy# Crear matriz simetrica
M = np.random.rand(4, 4)
A = (M + M.T) / 2

vals, vecs = np.linalg.eigh(A)
# vals: reales, ORDEN ASCENDENTE (diferente a eig)
# vecs: ortonormales: vecs.T @ vecs = I

# PCA manual con eigh:
X = np.random.rand(100, 5)
C = np.cov(X.T)                        # covarianza simetrica
vals, vecs = np.linalg.eigh(C)
# Reordenar descendente (eigh da ascendente)
idx = np.argsort(vals)[::-1]
components = vecs[:, idx]
variance_ratio = vals[idx] / vals.sum()

# Proyectar en k componentes principales
k = 2
X_pca = X @ components[:, :k]

copyA = np.array([[4,2],[1,3]])

# General (puede ser complejo)
np.linalg.eigvals(A)     # [5., 2.]

# Simetrica (siempre real, orden ascendente)
S = A + A.T
np.linalg.eigvalsh(S)

# Casos de uso comunes:

# Verificar definida positiva
es_spd = np.all(np.linalg.eigvalsh(S) > 0)

# Radio espectral (convergencia iterativa)
rho = np.max(np.abs(np.linalg.eigvals(A)))

# Numero de condicion espectral
eigs = np.linalg.eigvalsh(S)
cond = eigs.max() / eigs.min()

# Comprobar convergencia de metodo iterativo
converge = rho < 1.0

copyv = np.array([3.0, 4.0])

np.linalg.norm(v)              # L2 Euclidiana    = 5.0
np.linalg.norm(v, ord=1)      # L1 Manhattan     = 7.0
np.linalg.norm(v, ord=2)      # L2 (= default)   = 5.0
np.linalg.norm(v, ord=np.inf) # L-inf max(|v|)   = 4.0
np.linalg.norm(v, ord=0)      # L0 no-ceros      = 2.0
np.linalg.norm(v, ord=-np.inf)# min(|v|)         = 3.0

copyA = np.array([[1,2],[3,4]])

np.linalg.norm(A)              # Frobenius (default)
np.linalg.norm(A, ord='fro')  # Frobenius: sqrt(sum|aij|^2)
np.linalg.norm(A, ord=1)      # max suma columnas
np.linalg.norm(A, ord=np.inf) # max suma filas
np.linalg.norm(A, ord=2)      # valor singular maximo
np.linalg.norm(A, ord='nuc')  # norma nuclear: suma sigma_i

# Versiones modernas broadcast-aware
np.linalg.vector_norm(A, axis=1)    # norma L2 de cada fila
np.linalg.matrix_norm(A, ord='fro') # norma matricial

# Normalizar un vector
v_unit = v / np.linalg.norm(v)

# Normalizar filas de una matriz
A_norm = A / np.linalg.norm(A, axis=1, keepdims=True)

copyA = np.array([[1,2],[3,4]])

np.linalg.det(A)   # -2.0  (1*4 - 2*3)

# Para matrices grandes, det puede overflow/underflow
# slogdet: retorna (signo, log|det|) — siempre estable
sign, logdet = np.linalg.slogdet(A)
det = sign * np.exp(logdet)

# Batch: stacked matrices
batch = np.random.rand(100, 3, 3)
dets = np.linalg.det(batch)   # (100,)

# Digitos de precision perdidos por condicionamiento
import math
perdidos = math.log10(np.linalg.cond(A))

copyA = np.array([[1,2],[1.001,2]])   # casi singular

np.linalg.cond(A)            # norma 2 (default)
np.linalg.cond(A, p=1)       # norma 1
np.linalg.cond(A, p=np.inf)  # norma infinito
np.linalg.cond(A, p='fro')   # norma Frobenius

# Interpretacion:
# cond = 1        bien condicionada
# cond = 1e6      pierde ~6 digitos de precision
# cond > 1/eps    numericamente singular
# eps machine float64 = 1e-16

# Clasificacion practica
c = np.linalg.cond(A)
estado = "bien" if c < 1e4 else ("mal" if c < 1e12 else "singular")

copyA = np.array([[1,2,3],[4,5,6],[7,8,9]])
# fila3 = fila1 + fila2 -> rango 2

np.linalg.matrix_rank(A)              # 2
np.linalg.matrix_rank(A, tol=1e-10)   # tolerancia custom
np.linalg.matrix_rank(A, hermitian=True)  # mas rapido si simetrica

m, n = A.shape
rango_completo = np.linalg.matrix_rank(A) == min(m, n)
nulidad = n - np.linalg.matrix_rank(A)

copynp.trace(A)              # suma diagonal: 15
np.trace(A, offset=1)    # super-diagonal: 8
np.trace(A, offset=-1)   # sub-diagonal: 12
np.einsum('ii', A)       # equivalente con einsum

# Identidad: trace(A) = sum(eigvals(A))
# Norma Frobenius via traza:
frob = np.sqrt(np.trace(A.T @ A))

copyA = np.array([[2,1,-1],[-3,-1,2],[-2,1,2]], dtype=float)
b = np.array([8, -11, -3], dtype=float)

x = np.linalg.solve(A, b)    # [2., 3., -1.]
np.allclose(A @ x, b)         # True

# Multiples vectores b simultaneamente
B = np.random.rand(3, 5)     # 5 sistemas a la vez
X = np.linalg.solve(A, B)    # (3,5)

# Batch: stacked matrices
A_batch = np.random.rand(100, 4, 4)
b_batch = np.random.rand(100, 4)
x_batch = np.linalg.solve(A_batch, b_batch)  # (100,4)

copy# Lento e inestable
x = np.linalg.inv(A) @ b

# Correcto: 3x mas rapido y mas estable
x = np.linalg.solve(A, b)

# Multiples b con misma A: usar scipy LU
from scipy.linalg import lu_factor, lu_solve
lu, piv = lu_factor(A)
for b_i in many_bs:
    x_i = lu_solve((lu, piv), b_i)  # reutiliza LU

copy# Sistema sobredeterminado: mas ecuaciones que incognitas
A = np.array([[1,1],[1,2],[1,3],[1,4]])  # (4,2)
b = np.array([2, 3, 5, 4])              # (4,)

x, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None)
# x:         solucion LS
# residuals: ||Ax - b||^2 (vacio si rank < n)
# rank:      rango numerico de A
# sv:        valores singulares de A

# Regresion lineal clasica
t = np.linspace(0, 10, 100)
y = 3*t + 1 + np.random.randn(100)
A_reg = np.column_stack([t, np.ones_like(t)])
(slope, intercept), *_ = np.linalg.lstsq(A_reg, y, rcond=None)

copyA = np.array([[1,2],[3,4]])

# Inversa exacta (solo cuadradas no-singulares)
A_inv = np.linalg.inv(A)
np.allclose(A @ A_inv, np.eye(2))   # True

# Pseudoinversa Moore-Penrose (forma general)
B = np.random.rand(5, 3)     # (m,n) cualquier forma
B_pinv = np.linalg.pinv(B)   # (3,5)

# Propiedades de la pseudoinversa:
np.allclose(B @ B_pinv @ B, B)           # True
np.allclose(B_pinv @ B @ B_pinv, B_pinv) # True

# Solucion LS via pinv (equivalente a lstsq)
x = np.linalg.pinv(A) @ b

# Con tolerancia de singularidad personalizada
np.linalg.pinv(A, rcond=1e-10)

copy# Causas de LinAlgError:
# - inv(A) con A singular
# - solve(A,b) con A singular
# - cholesky(A) con A no definida positiva
# - eig/eigh sin convergencia

# Patron de manejo robusto
def safe_solve(A, b):
    try:
        cond = np.linalg.cond(A)
        if cond > 1e12:
            # Caer a minimos cuadrados con regularizacion
            return np.linalg.lstsq(A, b, rcond=None)[0]
        return np.linalg.solve(A, b)
    except np.linalg.LinAlgError:
        # Regularizacion Tikhonov como fallback
        n = A.shape[0]
        lam = 1e-6
        return np.linalg.solve(A.T @ A + lam*np.eye(n), A.T @ b)

# Verificar SPD antes de Cholesky
def safe_cholesky(A, jitter=1e-6):
    try:
        return np.linalg.cholesky(A)
    except np.linalg.LinAlgError:
        A_reg = A + jitter * np.eye(A.shape[0])
        return np.linalg.cholesky(A_reg)

copyA = np.array([[1,2,3],[4,5,6],[7,8,9]])

# Extraer diagonal
np.diagonal(A)              # [1,5,9]
np.linalg.diagonal(A)       # version broadcast-aware
np.diagonal(A, offset=1)    # [2,6] superdiagonal
np.diagonal(A, offset=-1)   # [4,8] subdiagonal

# Crear matrices diagonales
np.diag([1,2,3])            # matriz 3x3 diagonal
np.diag([1,2], k=1)         # superdiagonal
np.diag([1,2], k=-1)        # subdiagonal

# Transpuesta
A.T                          # transpuesta
A.conj().T                   # Hermitiana (conjugada transpuesta)
np.linalg.matrix_transpose(A) # version broadcast-aware

# Batch transpose: (N,m,n) -> (N,n,m)
batch = np.random.rand(10, 3, 4)
batch.transpose(0, 2, 1)     # (10,4,3)

copy# Identidad
np.eye(3)                    # I 3x3
np.eye(3, 4)                 # rectangular 3x4
np.eye(4, k=1)               # super-diagonal

# Ceros, unos
np.zeros((3,4))
np.ones((2,2))

# Simetrica definida positiva (SPD)
M = np.random.rand(4, 4)
SPD = M.T @ M + 1e-6 * np.eye(4)

# Matriz ortogonal (via QR)
M = np.random.randn(4, 4)
Q, _ = np.linalg.qr(M)     # Q es ortogonal
np.allclose(Q.T @ Q, np.eye(4))  # True

# Triangular inferior/superior
np.tril(A)                   # triangular inferior
np.triu(A)                   # triangular superior
np.tril(A, k=-1)             # estrictamente inferior